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The  International  Workshop  on  Computational  Electronics  was  held  in  Tempe, 
Arizona,  on  October  30  -  November  2,  1995.  This  workshop  was  the  fourth  in  a  series 
and  was  held  under  the  auspices  of  the  National  Center  for  Computational  Electronics  at 
the  University  of  Illinois  (K.  Hess,  Director).  Previous  conferences  have  been  held  in 
Urbana,  Leeds  (UK),  and  Portland.  In  general,  the  series  of  workshops  covers  all  aspects 
of  advanced  simulation  and  modeling  of  electronic  transport  in  semiconductors  and 
semiconductor  devices,  particularly  those  aspects  that  utilize  intensive  computational 
and/or  visualization  resources. 

The  workshop  is  intended  to  be  an  international  forum  for  the  discussion  of  the 
current  trends  and  the  future  directions  of  computational  electronics.  The  IWCE  attempts 
to  bring  together  engineers,  physicists,  computational  scientists,  and  applied 
mathematicians  in  the  area  of  simulation  and  modeling.  The  major  topics  of  this  past 
year’s  workshop  continued  to  be: 

•  Advances  in  two-  and  three-dimensional  simulations  utilizing  standard  techniques 
such  as  drift-diffusion  and  hydrodynamic  equations. 

•  Particle  simulation  methods  such  as  ensemble  Monte  Carlo,  molecular  dynamics,  and 
cellular  automata. 

•  Simulation  of  optical  processes,  optoelectronic  and  electro-optic  devices  which  e.g. 
need  the  incorporation  of  Maxwell’s  equations.\ 

•  Quantum  transport  and  quantum  devices. 

•  High  performance  computing  techniques  for  computational  electronics,  including 
parallelization,  vectorization,  and  improved  numerical  algorithms. 


1 


A  major  subsidiary  theme  for  this  past  year’s  workshop  was  process  modeling  and 
simulation.  To  accentuate  this  theme,  a  special  session  was  held  with  several  invited 
papers  addressing  this  topic.  This  focused  upon  both  feature-scale  and  reactor-scale 
process  modeling  and  simulation. 

Site 

The  IWCE  was  held  at  the  Radisson  Tempe  Mission  Palms  Hotel  in  downtown 
Tempe,  Arizona.  This  hotel  is  a  few  minutes  walk  from  Arizona  State  University,  and  is 
only  a  few  minutes  from  Phoenix’s  Sky  Harbor  International  Airport. 

Program/Organizing  Committee 

The  Program  and  Organizing  Committee  was  composed  of  internationally  known 
scientists  in  the  disciplinary  areas  covered  by  the  conference.  This  year’s  committee  is: 

•  David  K.  Ferry,  Arizona  State  University,  Chair 

•  Karl  Hess,  University  of  Illinois,  Director,  NCCE 

•  Chris  Ringhofer,  Arizona  State  University,  Local  Arrangements 

•  Carl  Gardner,  Arizona  State  University,  Publications 

•  Robert  Dutton,  Stanford  University 

•  Massimo  Fischetti,  IBM  Watson  Research  Center 

•  Stephen  Goodnick,  Oregon  State  University 

•  Chihiro  Hamaguchi,  Osaka  University 
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•  Joseph  Jerome,  Northwestern  University 

•  Peter  Markowich,  University  of  Berlin 

•  Hiroshi  Mizuta,  Hitachi  Central  Research  Lab. 

•  Massimo  Rudan,  Bologna  University 

•  R.  Kent  Smith,  A.T.&T.  Bell  Laboratories 

•  Chris  Snowden,  Leeds  University 

In  addition  to  the  Program/Organizing  Committee,  the  IWCE  benefits  from  the 
guidance  and  insight  of  an  international  advisory  committee,  which  is  composed  of: 

•  Herb  Bennett,  National  Institutes  of  Science  and  Technology 

•  Larry  Cooper,  Office  of  Naval  Research 

•  William  Coughran,  A.T.&T.  Bell  Laboratories 

•  Wolfgang  Fichtner,  ETH,  Zurich 

•  William  Frensley,  University  of  Texas  at  Dallas 

•  M.  Fukuma,  NEC  Corporation 

•  Harold  Grubin,  SRA,  Inc. 

•  Gerald  lafrate.  Army  Research  Office 

•  Thomas  Kerkhoven,  University  of  Illinois 

•  Steve  Laux,  IBM  Watson  Research  Center 

•  Thomas  C.  McGill,  California  Institute  of  Technology 

•  Mark  Pinto,  A.T.&T.  Bell  Laboratories 

•  A1  Tasch,  University  of  Texas 
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•  Peter  Vogl,  Walter  Schottky  Institut 

•  Naoki  Yokoyama 

Invited  Speakers 

Approximately  10  invited  speakers  were  selected  to  preface  the  major  topics 
described  earlier,  including  the  special  session  on  process  modeling  and  simulation. 
These  invited  speakers  will  be  drawn  from  internationally  recognized  experts  in  the  field, 
and  will  help  to  provide  an  introduction  to  the  contributed  papers  in  each  session. 

Student  Support 

One  of  the  most  important  aspects  of  the  IWCE  is  the  opportunity  for  graduate 
students  pursuing  research  in  this  area  to  discuss  their  work  with  internationally 
recognized  participants  in  the  workshop.  This  was  aided  by  approximately  12  student 
stipends  to  help  defray  the  cost  of  attendance  at  the  workshop  by  the  students.  The 
availability  of  this  support  was  advertised  in  the  second  announcement  and  call  for 
papers,  and  preference  was  given  to  students  who  presented  papers  at  the  workshop. 

Proceedings 

The  proceedings  of  the  4th  IWCE  will  be  published  as  a  regular  issue  of  VLSI 
Design.  This  will  be  an  issue  of  approximately  300  pages,  as  all  manuscripts  will  be 
contained  in  a  single  issue.  Papers  will  be  submitted  at  the  beginning  of  the  workshop 
with  regular  refereeing  occuring  during  the  workshop. 
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Conference  Results 


In  the  following  pages,  the  program  of  the  workshop  is  presented.  That  is 
followed  by  a  list  of  the  attendees  to  the  workshop. 

Next  Workshop 

The  Fifth  workshop  is  scheduled  for  May-June  1997,  and  will  be  held  at  the 
University  of  Notre  Dame,  Notre  Dame,  Indiana.  It  will  be  organized  by  Prof.  Wolfgang 
Porod. 
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Program 


MAI.  Quantum  Kinetic  Transport  under  High  Electric  Fields,  Nobuyuki  Sana  and  Akira  Yoshii, 
N7TLSI  Laboratories,  3-J  Morinisato-Wakamiwa,  Atsugi-shi,  Kanagawa  243’0I,  Japan,  As  the  device  size 
continues  to  shrink,  various  quantum  effects  are  expected  to  play  an  important  role  in  determining  the 
device  characteristics.  Among  others,  the  dynamical  quantum  effects  such  as  collisional  broadening  (CB) 
and  intracollisional  field  effect  (ICFE)  inevitably  come  into  play  in  conventional  ultrasmall  devices.[l]  A 
clear  understanding  of  these  effects  usually  requires  extensive  numerical  computations  and  is  an 
indispensable  matter  from  the  view  point  of  not  only  fundamental  physics  but  also  of  device  applications. 
In  the  present  talk,  we  will  give  a  brief  review  of  the  dynamical  quantum  effects  (CB  and  ICFE)  under  high 
electric  fields  and  show  our  recent  results  of  Monte  Carlo  calculations  of  electron  transport  for  3-,  2-,  and 
i-D  electron  gases.  We  derive  the  quantum  kinetic  (Barker-Feny)  transport  equation  from  the  reduced 
density  matrix  of  electrons  by  employing  the  projection  operator  approach.  It  is  shown  that  both  CB  and 
ICFE  are  closely  related  to  the  finite  collision  duration  time  with  phonons,  i.e.,  the  finite  life-time  of  the 
electron  yields  the  broadening  of  the  spectral  density  (CB)  and  the  energy  gain  (loss)  from  (to)  the  external 
electric  field  during  collision  duration  skews  as  well  as  broadens  the  spectral  density  (ICFE).  It  is  also 
demonstrated  that  ICFE  is  insignificant  in  bulk  up  to  extremely  large  electric  fields  (several  hundreds 
kV/cm),  whereas  it  becomes  more  and  more  important  even  under  moderate  electric  field  strengths  as  the 
dimensionality  of  the  structure  is  reduced.  This  is  because  ICFE  is  most  effective  when  the  momentum 
transfer  of  the  electron  due  to  the  scattering  with  phonons  is  directed  to  the  electric  field  and  the 
momentum  transfer  is  strictly  restricted  to  reduced  dimensions  in  lower  dimensional  structures.  More 
specifically,  in  the  case  of  a  1-D  rectangular  quantum  wire  in  GaAs  under  the  extreme  quantum  limit,  the 
total  electron-phonon  scattering  rate  becomes  asymmetric,  as  shown  in  Fig.  1,  with  respect  to  the  direction 
of  the  electron  propagation.  The  double  peak  structure  near  the  threshold  energy  (£  =  36  meV)  for  phonon 
emission  is  ascribed  to  the  shift  of  the  energy  detuning  due  to  ICFE.  Namely,  when  the  electron  travels  in 
the  opposite  (same)  direction  to  the  electric  field,  the  electron  gains  (loses)  energy  from  (to)  the  electric 
field  during  collision  duration  and  the  phonon  energy  is  effectively  reduced  (increased).  The  change  of  the 
effective  phonon  energy  is  about  4  meV  for  F  =  500  V/cm  and  even  greater  than  the  damping  factor  (~  2.5 
meV)  often  assumed  for  CB.[2]  Figure  2  shows  the  electric  field  dependence  of  the  drift  velocity  obtained 
from  the  Monte  Carlo  simulations  for  T  =  300  K.  It  should  be  noted  that  the  drift  velocity  is  greatly 
reduced  when  ICFE  is  taken  into  account  in  I-D  quantum  wires.  In  short,  extensive  numerical 
computations  are  successfully  employed  to  extract  a  significance  of  the  dynamical  quantum  effects  in  low 
dimensional  electron  gases.  Especially,  ICFE  would  become  crucial  in  the  analyses  of  electron  transport 
when  the  dimensionality  of  the  structure  comes  into  play. 

[1]  J.  R.  Barker  and  D.  K.  Ferry,  Phys.  Rev.  Lett.  42,  1779  (1979);  L.  Reggiani,  P.  Lugli,  and  A-P.  Jauho, 
J.  Appl.  Phys.  64,  3072  (1988). 

[2]  S.  Briggs,  B.  A.  Mason,  and  J.  P.  Leburton,  Phys.  Rev.  B  40,  12001  (1989);  J.  P.  Leburton  and  D. 
Jovanovic,  Semicond.  Sci.  Technol.  7,  B202  (1992). 


IV1A2.  A  Generalized  Tunneling  Formula  for  Quantum  Device  Modeling,  R.  Lake"*,  G.  Klimeck\  R,  C. 
Bowen,  D,  Jovanovic",  P.  Sotirelis",  and  W,  R,Frensley;  "Corporate  R&D,  Texas  Instruments 
Incorporated,  Dallas,  TX  75235,  ‘'Eric  Jonsson  School  of  Engineering,  University  of  Texas  at  Dallas, 
Richardson,  TX  75083-0688,  We  have  developed  novel  theory  and  numerical  algorithms  to  create  a 
general  purpose  quantum  device  simulator  which  includes  the  effects  of  charging,  scattering  (alloy 
disorder,  interface  roughness,  acoustic  phonon,  and  polar  optical  phonon),  and  a  number  of  bandstructure 
models  (2,  4,  8,  and  10  band  with  nearest  and  next  nearest  neighbor  coupling).  The  most  novel  and  useful 
invention  resulting  from  this  effort  is  a  generalized  tunneling  formula  based  on  generalized  boundary 
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conditions.  Tile  boundary  conditions  give  a  unified  treatment  of  emitter  quasi-bound  states  and  emitter 
continuum  states,  and  they  treat  thermalization  amd  scattering  induced  broadening  in  the  leads  with  no 
more  effort  than  the  standard  coherent  tunneling  approaches  [1].  The  single-band  version  of  the  theory  is 
briefly  described  in  references  [2,  3,  4]  and  single-based  numerical  results  based  on  the  theory  are 
described  in  references  [4,  5].  We  will  describe  both  the  single-band  and  multi-band  versions  of  the  theory 
and  we  will  give  both  single-band  and  multi-band  numerical  examples  of  device  simulations  using  the 
theory. 

[1]  D.  Landheer  and  G.  C.  Aers,  Superlattices  and  Microstructures  7,  17  (1990). 

[2]  R.  Lake,  in  Proceedings  of  the  Third  Internataonal  Workshop  on  Computational  Electronics  (Portland, 
OR,  May  18-20,  1994),  p.  239  . 

[3]  R.  Lake,  in  Quantum  Transport  in  Ultrasmall  Devices,  edited  by  D.  K.  Ferry,  H.  L.  Grubin,  and  C. 
Jacoboni  (Plenum  Press,  New  York,  1995),  pp.  521-524. 

[4]  G.  Klimeck,  R..  Lake,  R.  C.  Bowen,  W.  R.  Frensley,  and  T.  Moise  Quantum  Device  Simulation  with  a 
Generalized  Tunneling  Formula,  (preprint,  April  1995) . 

[5]  G.  Klimeck,  R.  Lake,  R.  C.  Bowen,  W.  R.  Frensley,  and  D.  Blanks,  in  the  1995  53rd  Annual  Device 
Research  Conference  Digest,  p.  52,  Charlottesville,  VA,  June  19-21, 1995. 


MA3.  Monte  Carlo  Simulation  of  HEMT  Using  Self-Consistent  Method,  K  Ueno,  S.  Yamakawa  and  C. 
Hamaguchi,  Dept  of  Electronic  Engineeringy  Faculty  of  Engineering,  Osaka  University,  Suita  City,  Osaka 
565,  Japan;  and  K.  Miyatsuji,  Electronics  Research  Laboratory,  Matsushita  Electronics  Corporation, 
Takatsuki  City,  Osaka  569,  Japan,:  We  report  here  Monte  Carlo  simulation  of  HEMT  by  using  self- 
consistent  method,  where  the  real  device  structures  are  taken  into  account.  The  electronic  states  of  2DEG  in 
the  channel  are  calculated  by  solving  Schrodinger  equation  and  Poisson  equation  self-consistently,  by 
taking  the  F,  L  and  X  valleys,  where  the  electrons  in  the  three  valleys  are  assumed  to  be  two-dimensional. 
The  eigen  states  of  2DEG  in  the,  F  L  and  X  valleys  obtained  by  self-consistent  calculations  are  used  to 
evaluate  the  scattering  rates.  For  the  scattering  of  the  electrons,  we  take  account  of  the  acoustic  phonon, 
polar  optical  phonon,  intervalley  phonon  and  ionized  impurity  scatterings,  where  the  nonparabolic  energy 
structures  are  used.  In  our  previous  simulation  of  HEMT  the  electric  field  parallel  to  channel  is  treated  as 
uniform,  resulting  in  unnegligible  error  in  drain  current  of  saturation  region.  In  the  real  device  the  sheet 
electron  density  in  the  channel  is  not  homogeneous  and  the  electric  field  along  the  channel  is  not  uniform, 
higher  near  the  drain.  In  the  present  simulation  we  take  into  account  of  the  nonuniformity  of  the  electric 
field  due  to  the  sheet  electron  density  distribution.  The  electric  field  at  each  point  of  a  real  device  is 
evaluated  by  calculating  two-dimensional  Poisson  equation  where  the  distribution  of  the  electrons  and  the 
ionized  impurities  are  taken  into  account.  At  high  electric  fields  electrons  are  heated  to  form  hot  electrons 
which  behave  as  three-dimensional  electrons  rather  than  2DEG.  In  addition  hot  electrons  make  real  space 
transfer  between  the  GaAs  layer  and  the  AlGaAs  barrier.  These  effects  are  taken  into  account  by  assuming 
that  an  electron  with  energy  higher  than  the  barrier  after  the  scattering  is  three-dimensional.  Such  a  three- 
dimensional  electron  is  staying  in  the  GaAs  layer  or  is  transferred  into  the  AlGaAs  layer  depending  on  the 
wave  vector  of  the  electron.  Simulation  were  carried  out  for  a  typical  HEMT  structure  with  sheet  electron 
density  5  x  10“cm"^  and  A1  content  x  =  0.3.  Obtained  results  at  300K  will  be  reported. 


MA4.  The  Smooth  Effective  Potential  for  the  Quantum  Hydrodynamic  Model,  Carl  L  Gardner  and 
Christian  Ringhofer,  Department  of  Mathematics, Arizona  State  University,  Tempe,  AZ  85287- J 804.  An 
"extended"  quantum  hydrodynamic  (QHD)  model,  valid  for  potentials  with  discontinuities,  is  derived  in  the 
Bom  approximation  to  the  Bloch  equation.  The  quantum  potential  appearing  in  the  stress  tensor  Pij  and 
energy  density  W  is  valid  to  all  orders  of  ^  ^  and  to  first  order  in  f5  V,  and  involves  both  a  smoothing 
integration  of  the  classical  potential  over  space  and  an  averaging  integration  over  temperature.  At  potential 
barriers  in  semiconductors,  the  new  effective  stress  tensor  and  energy  density  are  more  tractable  analytically 
and  numerically  than  in  the  original  O(h^)  QHD  theory.  By  cancelling  the  leading  singularity  in  the 
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classical  potential  at  a  barrier  and  leaving  a  residual  smooth  effective  potential  (with  a  lower  potential 
height)  in  the  barrier  region,  the  effective  stress  tensor  makes  the  barrier  partially  transparent  to  the  particle 
flow  and  provides  the  mechanism  for  particle  tunneling  in  the  QHD  model.  Numerical  comparisons 
demonstrate  excellent  agreement  between  the  first  three  moments  of  the  equilibrium  full  density  matrix  and 
the  effective  density  matrix  for  the  Bloch  equation  with  a  barrier  potential  for  fiAV  <  1,  and  good 
qualitative  agreement  for  jSAV  <  10.  The  moments  of  the  0{tl  density  matrix  are  in  severe  quantitative 
and  qualitative  disagreement  with  the  moments  of  the  full  density  matrix.  Mathematical  results  on  the 
convergence  of  our  iteration  method  for  solving  the  Bloch  equation  and  numerical  simulations  of  a  resonant 
tunneling  diode  comparing  solutions  of  the  extended  and  0(^  QHD  equations  are  also  presented. 


MAS.  Quantum  Transport  Simulation  of  the  DOS  Function,  Self-Consistent  Fields  and  Mobility  in 
MOS  Inversion  Layers*,  Dragica  Vasileska,  Terry  Eldndeg@,  Paolo  Bordone^,  and  David  K.  Ferry, 
Center  for  Solid  State  Electronics  Research,  Arizona  State  University,  Ternpe,Az,  85287-6206.  We 
propose  a  simulation  of  the  self-consistent  fields  and  mobility  in  (1(X))  Si-inversion  layers  for  arbitrary 
inversion  charge  densities  and  temperatures.  Green's  functions  are  employed  in  obtaining  the  analytic 
expressions  for  the  broadening  of  the  states,  real  shift  in  the  subband  energies  and  conductivity.  The  initial 
potential  energy  profile  is  calculated  numerically  by  a  variational  approach.  The  Poisson,  Schrodinger  and 
Dyson  equation  for  the  retarded  Green’s  functions  are  solved  independently  for  the  corresponding 
unknowns  and  the  self-consistency  is  achieved  through  the  outer  iterative  loop,  split  up  into  two  parts. 
First,  the  Schrodinger  and  Poisson  equations  are  solved  by  assuming  an  ideal  density  of  states  (DOS) 
function.  In  the  second  loop,  we  use  the  real  DOS  function.  The  potential  energy  profile  for  the  next 
iteration  is  calculated  by  using  a  combination  of  fixed-  and  extrapolated-convergence  factor  scheme.  The 
exchangecorrelation  corrections  are  calculated  by  using  the  density-functional  formalism.  The 
wavevectordependent  matrix  elements  for  the  various  scattering  mechanisms  (surface-roughness  scattering, 
depletion  layer  and  interface/oxide  Coulomb  scattering,  acoustic  and  non-polar  optical  phonon  scattering) 
are  calculated  within  the  self-consistent  Bom  approximation.  Screening  is  treated  within  the  RPA.  The 
diagonal  polanzation  function  is  used  in  the  calculation  of  the  screened  matrix  elements.  The  simulation 
results  suggest  that  the  proposed  theoretical  model  gives  mobilities  which  are  in  very  good  agreement  with 
the  experimental  data. 

*  Work  supported  by  the  Office  of  Naval  Research 

@  McDonnel  Douglass  Helicopter  Systems,  5000  E.  McDowell  Rd,  Mesa,  Az,  85215. 

#  Permanent  address:  Dipartimento  di  Fisica  ed  Institute  Nazionale  di  Fisica  della  Materia,  Universita  di 
Modena,  Via  Campi  213/A,  41  IflO  Modena,  Italy. 

[1]  D.K.  Ferry,  Phys.  Rev.  B  14, 5364  (1976). 

[2]  S.  Takagi,  M.  Iwase,  and  A.  Toriumi,  lEDM  Tech,  Dig.  398  (1988). 


MA6.  Study  of  Interface  Roughness  Dependence  of  Electron  Mobility  in  Si  Inversion  Layer  using 
Monte  Carlo  Method,  S.  Yamakawa,  H.  Ueno  and  C.  Hamaguchi,  Department  of  Electronic  Engineering, 
Faculty  of  Engineering,  Osaka  University,  Suita  City,  Osaka  565,  Japan;  K.  Masaki,  Anan  College  of 
Technology,  Anan  City,  Tokashima  744,  Japan;  K.  Miyatsuji,  Electronics  Research  Laboratory,  Matsushita 
Electronics  Corporation,  Takatsuki  City,  Osaka  569,  Japan;  U.  Ravaioli,  Department  of  Electrical  and 
Computer  Engneering,  University  of  Illinois  at  Urbana-Champaign,  Urbana,  IL  61801.  We  report  the 
effect  of  interface  roughness  scattering  on  the  electron  mobility  in  MOSFET.  Mobility  of  the  two- 
dimensional  electron  gas  (2DEG)  in  MOSFETs  is  one  of  the  most  important  parameters  from  the  viewpoint 
of  device  operation.  Since  a  strong  gate  field  results  in  confinement  of  electrons  at  the  si/si02  interface  in  a 
MOSFET,  with  the  formation  of  a  2DEG,  the  electron  transport  in  the  inversion  layer  is  known  to  be 
different  from  the  case  of  threedimensional  carriers  in  bulk  materials.  So  first  we  calculate  the  subband 
stmctures  by  solving  Schrodinger  equations  and  Poisson  equations  self-consistently.  Since  at  strong  normal 
fields  the  electrons  are  confined  in  a  narrow  region  of  the  conduction  channel  near  the  interface,  there  is 
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considerable  scattering  due  to  roughness  of  the  Si/Si02  interface,  resulting  in  a  reduction  of  the  effective 
mobility.  Electron  mobility  in  the  inversion  layer  of  a  MOSFET,  formed  on  the  (100)  silicon  surface,  is 
calculated  by  using  a  Monte  Carlo  approach  which  takes  into  account  size  quantization,  acoustic  phonon 
scattering,  intervalley  phonon  scattering  and  surface  roughness  scattering.  Degeneracy  is  also  considered 
because  it  is  important  at  higher  normal  effective  fields  (high  gate  voltages).  The  main  emphasis  is  placed 
on  the  influence  of  the  specific  autocovariance  function,  used  to  express  the  surface  roughness,  on  the 
electron  mobility.  Here  we  compare  the  mobilities  obtained  using  exponential  and  Gaussian  autocovariance 
functions.  It  is  found  that  the  electron  mobility  calculated  with  roughness  scattering  rates  based  on  the 
exponential  function  shows  a  good  agreement  with  experiments.  The  effect  of  screening  on  the  roughness 
scattering  is  also  discussed. 


MA7.  Lattice  Effects  in  the  Complex  Subband  Dispersion  of  2DEG  Semiconductor  Waveguide 
Structures  Subject  to  a  Perpendicular  Magnetic  Field*,  G.  Edwards  and  D.  K.  Ferry,  Center  for  Solid 
State  Electronics  Research,  Arizona  State  University,  Tempe,  AZ,  85287-6206.  In  modeling  2DEG 
magneto-transport  experiments  it  is  important  to  have  knowledge  of  the  electronic  states  subject  to  a 
magnetic  field  perpendicular  to  the  plane  of  the  2DEG,  and  the  waveguide  confinement  potential.  The 
Schrodinger  equation  including  the  B  field  and  the  confinement  potential  can  be  solved  by  a  discretization 
procedure  and  hence  putting  the  wavefunction  'field'  on  a  lattice.  When  the  lattice  constant  is  much  smaller 
than  the  Fermi  wavelength  the  lattice  model  should  be  able  to  reproduce  the  true  continuum  situation,  for 
states  up  to  the  Fermi  level.  We  present  numerical  results,  within  the  lattice  model,  for  the  full  edge  state 
(magneto-electric  states)  complex  subband  dispersion  of  a  rectangular  waveguide,  including  both  the  real 
'bands'  and  the  complex  evanescent  'bands’.  The  full  complex  dispersion  is  needed  in  treating  a 
heterogeneous  structure  such  as  a  ballistic  cavity  or  disordered  quantum  wire,  in  a  B  field,  when  a 
rectangular  waveguide  section  is  used  as  a  lead  region  to  inject  current.  The  waveguide  states  subject  to  the 
B  field  can  be  obtained  analytically  by  perturbation  theory  or  WKB  theory.  At  a  fairly  fine  degree  of 
discretization  the  form  of  our  numerical  purely  real  subband  solutions  'agree'  well  with  the  analytic  real 
band  solutions.  However,  far  from  the  bandedges,  our  numerical  evanescent  solutions  can  have  a  different 
topology  to  the  analytic  evanescent  solutions.  We  find  that  a  very  fine  level  of  discretization  is  necessary  to 
describe  the  evanescent  states  accurately,  deep  into  the  complex  dispersion  region  . 


*  Work  supported  by  NEDO  under  the  International  Joint  Research  Program  'Nanostructures  and  Electron 
Waves  Project'. 


PI,  Simulation  of  a  Single  Electron  Tunnel  Transistor  with  Inclusion  of  Inelastic  Macroscopic 
Quantum  Tunneling  of  Charge,  Christoph  Wasshuber  and  Hans  Kosina,  Institute  for  Microelectronics, 
TU-Vienna,  Gusshausstrasse  27-29/E360,  A-1040  Wien,  Austria.  Until  now  Single  Electron  Tunnel  (SET) 
devices  were  simulated  by  either  neglecting  macroscopic  quantum  tunneling  of  charge  (q-MQT)  or 
approximating  it.  Thus,  we  simulated  a  SET  transistor  with  the  full  non-approximative  inclusion  of  inelastic 
q-MQT  or  cotunneling.  A  Monte  Carlo  method  was  used  to  simulate  electrons  that  tunnel  back  and  forth 
through  the  two  tunnel  junctions  of  the  SET  transistor  and  co-tunnel  back  and  forth  through  both  junctions 
simultaneously.  In  the  coulomb  blockade  regime  and  at  low  temperature  the  q-MQT  effect  dominates  the 
current  through  the  transistor.  The  thermally  agitated  normal  tunneling  is  orders  of  magnitude  smaller. 
Resonances  in  the  I-V  characteristic  were  found.  The  resonant  peaks  decrease  with  increasing  temperature. 
This  resonance  does  not  originate  from  the  normal  tunnel  effect,  like  in  the  well  known  resonant  tunneling 
in  double  barriers,  but  from  the  co-tunnel  effect.  Thus  the  simulation  shows  a  new  feature  of  the  SET 
transistor,  that  can  be  helpful  for  particular  measurements  of,  for  example,  the  coulomb  energy  or  related 
capacitances,  or  can  be  exploited  in  new  devices.  This  resonance  is  not  yet  experimentally  verified. 


P2.  Wireless  Single-Electron  Logic  Biased  by  Alternating  Electric  Field,  Alexarlder  N.  Korotkov, 
Department  of  Physics,  State  University  of  New  York  at  Stony  Brook,  NY  11794-3800.  Single-electron 
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effects  in  systems  of  small-capacitance  tunnel  junctions  can  possibly  be  used  as  the  physical  basis  for  a 
new  generation  of  ultradense  digital  electronics.[l]  Logic/memory  based  on  single-electron  transistors[l,2] 
and  logic/memory  which  uses  single  electrons  to  represent  digital  information[l]  have  been  discussed  in 
the  literature.  In  the  present  work  a  new  type  of  single-electron  logic  is  proposed. [3]  The  basic  element  is 
a  short  chain  of  islands  which  shows  bistable  polarization  and  affects  the  polarization  of  neighboring 
chains.  In  contrast  to  previous  approaches  wires  are  not  used  (which  is  a  very  favorable  feature  when 
considering  nanometer-size  circuits  and  possibly  molecular  structures),  and  the  circuits  are  biased  by  an 
external  electric  field.  The  proposed  logic  has  some  similarity  to  circuits  suggested  for  "ground  state 
computing".[4]  The  principal  difference  is  that  we  use  a  conventional  dissipative  method  of  switching 
elements  (power  is  supplied  by  an  alternating  electric  field)  which  makes  the  logic  simple,  robust,  and 
relatively  fast  and  allows  the  circuits  consisting  of  arbitrary  large  number  of  gates.  Monte-Carlo  simulation 
of  logic  gates,  propagation  lines,  and  fan-out  circuits  for  a  particular  geometry  shows  parameter  margins  of 
about  5%.  This  number  depends  on  the  chosen  layout  (positions  and  sizes  of  conducting  islands)  and  can 
probably  be  increased  up  to  10-15%. 

The  work  was  supported  by  AFOSR  and  ONR/ARPA. 

[1] .  D.  V.  Averin  and  K.  K.  Likharev,  in:  Single  Charge  Tunneling,  ed.  by  H.  Grabert  and  M.  Devoret 
(Plenum,  New  York,  1992),  p.  311. 

[2] .  A.  N.  Korotkov,  R.  H.  Chen,  and  K.  K.  Likharev,  J.  Appl.  Phys.  78  (August  15,  1995). 

[3]  A.  N.  Korotkov,  to  be  published. 

[4] .  P.  D.  Tougaw,  C.  S.  Lent,  and  W.  Porod,  J.  Appl.  Phys.  74,  3558  (1993). 


P3.  Single-Electron  Parametron,  Konstantin  K.  Likharev  and  Alexander  N.  Korotkov, 

Department  of  Physics,  State  University  of  New  York  at  Stony  Brook,  NY  11794-3800,  Since  1987,  several 
families  of  single-electron  logic  devices  have  been  suggested[l].  In  these  devices,  digital  bits  are  presented 
by  single  electrons,  which  are  manipulated  using  the  Coulomb  blockade  effect.  Most  of  the  devices 
suggested  earlier  imply  the  use  of  long  wires  for  delivery  of  the  power  supply  to  the  logic  gates.  In  this 
work  we  propose  a  new  type  of  single-electron  logic  consisting  of  three-island  elements  biased  by  a  rotating 
electric  field  which  also  plays  the  role  of  a  global  clock.  At  a  certain  phase  of  field  rotation  the  cell 
becomes  polarized  by  the  transfer  of  a  single  electron  from  the  central  island  to  one  of  the  edge  islands. 
The  symmetry  of  the  situation  is  broken  by  the  electric  field  created  by  earlier  polarized  neighboring  cells. 
The  idea  is  very  similar  to  that  used  in  ac  and  dc  parametrons  (see,  e.g..  Ref.  3  and  references  therein),  so 
we  call  our  basic  cell  the  single-electron  parametron.  In  contrast  to  the  recently  suggested  "wireless  single¬ 
electron  logic”[2],  the  parametron  is  physically  reversible.  As  a  consequence,  for  example,  the  shift  register 
based  on  these  cells  can  dissipate  power  less  than  ksT  per  element  per  switching.  We  have  carried  out 
extensive  analytical  calculations  and  numerical  simulations  of  the  single-electron  parametron  on  both 
equivalent-circuit  and  geometrical-model  levels.  Results  for  parameter  margins,  power  consumption  and 
operation  speed  will  be  presented  at  the  meeting. 


The  work  was  supported  by  AFOSR  and  ONR/ARPA. 

[1] .  D.  V.  Averin  and  K.  K.  Likharev,  in:  Single  Charge  Tunneling,  ed.  by  H.  Grabert  and  M.  Devoret 
(Plenum,  New  York,  1992),  p.  311. 

[2]  A.  N.  Korotkov,  to  be  published. 

[3] .  K.  K.  Likharev,  IEEE  Trans.  Magn.  13, 245  (1977). 


P4.  Parallel  Computation  for  Electronic  Waves  in  Quantum  Corrals*  Henry  K,  Harburyt  and 
Wolfgang  Porod^,  University  of  Notre  Dame,  ‘Science  Computing  Facilities,  *Dept,  of  Electrical 
Engineering,  Notre  Dame,  IN  46556,  Recent  scanning  tunneling  microscopy  (STM)  studies  on  the  (111) 
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faces  of  noble  metals  have  directly  imaged  electronic  surface-confined  states.  Dramatic  standing-wave 

on  Cu{ni).[l,2]  We  solve  for  the  local  density  of  electronic  states  (LDOS)  in  these  structures  using  a 
coherent  elastic  scattering  theory.  We  seek  solutions  of  the  two-dimensional  Schrodinger  equation 
compatible  with  Sommerfeld  radiation  conditions[3.4]  in  these  "leaky"  confinement  structures^  We  Lsent 

whll  calculations  to  the  reported  experimental  differential  conductance  STM  data 

hich  revea  exce  ent  agreement  with  our  elastic  scattering  theory.  We  also  investigate  possible  "surface 
waveguides",  similar  to  split-gate  type  structures,  by  modeling  a  wide-narrow-wl  Ta adatm 
confineinent  structure  for  the  Cu(lll)  surface.  We  implement  non-reflecting  boundary  conditions  thTch 
ymptoucally  satisfy  the  Sommerfeld  radiation  condition  on  the  boundary  using  a  local  method  based  on 
the  work  of  Bayhss  and  coworkers.[3.4]  TTie  large  matrices  generated  by  the'discrlTtlfn  oS^ 
quantum  corral  structures  require  the  use  of  sparse  matrix  methods.  In  addition,  a  parallel  finite  element 
u  ion  was  undemken  using  the  message  passing  interface  standard  (MPI)  and  the  Portable  Extensible 


Th«  wk  has  bean  supponri  by  ihe  Office  of  Naval  Reseaseh  and  U,.  Air  Force  Office  of  Seienrifie 

[1]  M.  F.Crommie.  C.  P.  Lutz,  and  D.  M.  Eigler.  Science  262.  218  (1993) 

C.  P.  Lutz,  and  D.  M.  Eigler,  Nature  369,' 464  (1994) 

Portfan^'^Uiwi^  ’  Workshop  on  Computational  Electronics, 

[4]  A.  Bayliss,  M.  Gunzburger,  and  E.  Turkel,  SIAM  J.  Appl.  Math.  42  430  (1982). 


f  Electron  Transport  in  SUicon  nMOSFET  Inversion  Layers, W. -AT  Shih  S 
rt  J7’  •  and  A.F.  Tasch,  Jr.,  Microelectronic' Research  ^ 

■  t  A  ^  at  Austin,  Austin,  TX  78712.  Despite  decades  of  research  effort  [11  problems 
socmted  with  earner  transport  in  the  inversion  layers  of  silicon  MOSFETs  continue  to  capture  the 
»enb™  of  are  res^cb  coinmoni.y.r2-IJ  Today’s  advanced  devices  introduce  further  complicafi^ns  wfih 
oversS  Of  channel  length  that  introduces  pronounced  non-local  effects  such  as  velocity 

i  h  •  •"'^estigate  and  anticipate  future  device  performance  issues,  simulation  tools  thS 

rely  on  physically  rather  than  empincally  based  models  are  required.  We  present  here  a  single  nanieie 
Monte  Carlo  (MC)  simulation  on  uniform  silicon  nMOSFET  inversion  layers  within  and  beyond^th’e^ohmic 
regime  for  vanous  substrate  doping  levels  and  gate  biases.  Inversion  layer  electrons,  whS  mlroHir 

in  silicon,  form  a  two-dimensional  electron  gas 
Lica?‘of^ir"?^°”  described  by  a  two-dimensional  multi-subband  BTE  for  the  fiefds 

typ  al  of  the  linear  re^on.  We  obtain  the  subband  structure  of  the  2DEG  by  self-consistentlv  solving  the 
e  ective-mass  Scluodinger  equation  and  the  Poisson  equation  within  UTMINIMOS  [51.  an  advanced 
evice  simulator  that  encompasses  drift  diffusion,  hydrodynamic  and  MC  models  The  non 

^cture,  using  a  formalism  that  builds  on  the  work  ol  Fischetti  et  al.[3]  Scattering  in  the  2D  svstem  due  to 
phonons  ^d  surface  roughness  is  considered  using  the  formalism  dW"b^  W  e 
espectively  with  sekction  rules  for  the  intervalley  phonon  processes  established  by  Ferry  m  DegenerlS 
fieH  h  .The  universal  relationship  between  the  effective  mobility  and  the  ^sverse  Effective 

mobilkv  ^Srt&ment  between  the  simulated  and  experimental  low-field 

mobility  has  been  achieved  for  both  300K  and  77K.  The  simulated  results  are  compared  with  the 
expenmental  work  of  Takagi  [8]  at  77K  and  with  the  UT  mobility  model  [21  at  300K  T^rrrr 
model  h,^  been  ebamemriz«l  wid.  .  ,„ge  eumbe,  of  devLTv ^1“ 

.e  evemse  mobifi,  of  elecboo.  J  die 
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degenerate)  subbands  are  illustrated  in  Fig.  3.  At  high  lateral  field,  the  electron  drift  velocity  and  energy 
distribution  are  illustrated  in  Figures  4  through  6.  The  calculated  inversion  layer  saturation  velocity  is 
found  to  be  independent  of  the  transverse  field  and  agrees  well  with  experimental  data  as  described  by  the 
UT  mobility  model.  With  different  numbers  of  subbands  included  in  the  simulation,  the  extracted  electron 
energy  distribution  indicates  the  necessity  of  including  either  more  subbands  or  a  classical  simulation 
domain  if  the  tail  of  the  energy  distribution  is  desired.  In  summary,  we  have  developed  a  Monte  Carlo 
simulation  tool  that  accurately  reproduces  and  explains  the  experimental  universal  mobility  curve  and  high- 
field  transport  characteristics  in  uniform  silicon  nMOSFET  inversion  layers  at  different  temperatures.  The 
MC  tool  also  provides  access  to  important  but  experimentally  unavailable  microscopic  information. 

This  work  was  supported,  in  part,  by  the  Semiconductor  Research  Corporation  (SRC),  the  Texas  Advanced 
Technology  Program,  Motorola  and  AMD. 
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[3]  M.  V.  Fischetti  et  al.,  Phys.  Rev.  B,  vol.  48, 4,  p.  2244,  1993. 

[4]  F.  Assaderaghi  et  al.,  lEDM  Tech  Dig.,  p.  479,  1994. 

[5]  C.-F.  Yeap  et  al.,  NUPADS  V,  p.  15-18,  1994.  X.  L.  Wang  et.  al.,  J.  Appl.  Phys.,  Vol.  73,  7,  p.3339, 
1993. 

[6]  P.  J.  Price.  Anls.  Phys.  vol.  133,  p.  217,  1981.  Y.  C.  Cheng,  Surf.  Sci.  vol.  27,  p.  663,  1971. 

[7]  D  K.  Fcity,  Surf.Sci.,.  vol.  57,  p.  218, 1976. 
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P6.  A  Numerical  Study  of  Rare  Events  in  Single  Electron  Devices.  L  R.  C.  Fonseca,  A.  N.  Korotkov, 
and  K.  K.  Likharev;  Department  of  Physics,  State  Uniuersity  ofNeun  York  at  Stony  Brook,  NY  11794-3800. 
During  the  past  decade  several  devices  based  on  correlated  tunneling[l]  of  single  electrons  have  been 
suggested  and  tested.[2,3]  These  devices  may  be  used  in  various  analog  and  digital  systems,  in  particular 
as  ultradense  logic/memory  units[l]  and  as  fundamental  standards  of  dc  current.[4]  However,  tools  for 
their  computer-aided  design  which  would  account  for  all  the  most  important  physical  phenomena  have  not 
been  available  until  very  recently.  In  my  talk  I  will  describe  a  new  algorithm[5]  suitable  for  studies  of  the 
dynamics  of  single  electron  devices  presenting  arbitrary  combinations  of  small  tunnel  junctions, 
capacitances,  and  voltage  sources.  In  contrast  to  Monte-Carlo  schemes,  the  method  is  useful  both  to 
analyse  the  "classical"  behavior  of  the  system  and  for  calculating  small  deviations  from  this  behavior,  due 
to  the  finite  speed  of  applied  signals,  thermal  activation,  and  cotunneling.  I  will  show  results  from  the 
analysis  of  errors  in  memory  cells  and  current  deviations  in  proposed  dc  current  standards. [4]  I  will  also 
discuss  the  effect  of  background  charge  on  the  behavior  of  these  devices. 


This  work  was  supported  by  AFOSR,  and  CNPq-Brazil. 

[1]  D.V.  Averin  and  K.K.  Likharev  in,Single  Charge  Tunnling,.  ed.  by  H.  Grabert  and  M.  Devoret 
(Plenum,  New  York,  1992),  p.  311. 

[2]  J.M.  Martinis  et  al.,  Phys.  Rev.  Lett.  72,  904  (1994). 

[3]  P.D.  Dresselhaus  et  al.,  Phys.  Rev.  Lett.  72,  3226  (1994). 

[4]  H.  Pothier  et  al.,  Europhys.  Lett.  17, 249  (1992). 

[5]  L.R.C.  Fonseca  et  al.,  J.  Appl.  Phys.  79  (September  1,  1995). 


P7.  New  Method  for  Computing  Resonances  and  Transmission  Properties  for  Quantum  Devices,  V. 
Mandelshtam,  T.  Ravuri  and  H.  S.  Taylor;  Department  of  Chemistry,  University  of  Southern  California, 
Los  Angeles,  CA  90089-0482.  Two  relatively  new  methods,  the  Spectral  Projection  method  and  the 
Stabilization  Method,  of  implementing  scattering  calculations  are  described,  and  are  here  applied  to  two 
devices.  Both  methods  use  essentially  short  range  spectral  projection  operators  to  produce  a  complete  set 
of  solutions  of  the  wave  equation  that  need  be  valid  only  inside  the  interaction  region.  While  the  Spectral 
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Projection  method  is  more  generic  than  the  Stabilization  Method  which  is  based  on  using  the  more  difficult 
to  compute  spectral  density  operator,  the  latter  becomes  very  efficient  when  narrow  resonances  exist.  For 
problems  of  small  size  methods  are  practical  in  the  sense  that  they  involve  only  real,  symmetric  matrices 
resulting  from  Hamiltonians  represented  on  basis  sets.  For  more  challenging  larger  systems  the  Spectral 
Projection  method  lends  itself  to  a  very  efficient  time  independent  iterative  procedure  that  obtains 
transmission  results  simultaneously  at  all  energies.  This  procedure  uses  modified  Chebyshev  recursion 

relations  to  essentially  expand  the  operator  -  H  +  f  ,  where  f  is  an  absorbing  optical  potential  that 

eliminates  the  need  for  large  grids.  This  latter  procedure  requires  minimal  storage  and  the  resulting  series 
converges  rapidly  in  a  manner  that  is  uniform  in  energy.  For  small  systems  resonances  can  be  computed 
directly  using  stabilization.  "Directly"  means  no  S  matrix  or  transmission  coefficient  is  needed.  For  large 
systems  the  modified  Chebyshev  version  of  the  spectral  projection  operator  can  be  used  to  create  a  basis 

with  which  H  +  Fean  be  represented  and  diagonalized  to  yield  a  direct  calculation  of  resonance  energies. 
Large  problems  requiring  up  to  several  hundred  thousand  functions  or  grid  points  can  be  treated  using  this 
method. 


P8.  Quantum  Contributions  and  Violations  of  the  Classical  Law  of  Mass  Action,  H  L  Grubin  and  T. 
R.  Govindan;  Scientific  Research  Associates,  Inc.,  Glastonbury,  Connecticut..  The  law  of  mass  action 
refers  to  the  product  of  the  concentration  of  reactants,  where  for  a  particular  material  the  product  is  a 
function  of  temperature  only.  In  the  case  of  nondegenerate  semiconductors  this  is  represented  by  the 
product  of  electrons  and  holes: 

NP=A^.  =Nfi/^exp-Eg  AgT.  For  degenerate  material  NP  is  a  function  of  concentration  as  well  as 

temperature.  Since  the  work  of  Wigner  it  has  been  known  that  when  there  are  strong  concentration 
gradients  the  carrier  density  is  no  longer  given  by  a  simple  Boltzmann  like  distribution.  As  discussed  by 
Ancona  and  laffate,  for  non-degenerate  semiconductors  the  density  may  be  approximated  by  N  = 

exp'[Ep(j:)  +  (ft^/6m)V^VN/VNj/%jr  where  we  have  incorporated  the  quantum  potential  in  the 


exponential.  The  quantum  modification  of  the  density  suggests  that  the  product  of  NP  under 
thermodynamic  equilibrium  would  depend  on  features  other  than  temperature.  Certainly  this  would  be 
expected  at  heterostructure  interfaces,  even  if  effective  mass  variations  are  ignored.  The  question 
addressed  in  this  paper  is:  are  significant  variations  in  the  product  NP  to  be  expected  in  ostensibly 
classical  structures  such  as  a  bipolar  transistor.  To  answer  this  question  we  have  solved  the  quantum 

Liouville  equation  (avoiding  the  use  of  the  quantum  potential):  ih — - —  =  j^(f),p  (f)J  under 


conditions  of  thermodynamic  equilibrium.  The  equation  was  coupled  to  Poisson's  equation  and  solved  in 
the  coordinate  representation  for  two  generically  nondegenerate  classical  structures:  NIN  and  NPN 
devices,  where  in  all  cases  the  I  and  P  regions  were  each  50  nm  in  length.  Devices  of  200  nm  and  400  nm 
were  considered.  Several  calculations  were  performed  for  parameters  appropriate  to  GaAs;  parametric 
studies  were  also  performed.  In  these  studies  the  N  region  was  doped  to  10“/m’,  (which  is  too  high  to 
strictly  satisfy  the  nondegeneracy  condition).  Calculations  at  lower  densities  will  be  discussed  at  the 
meeting.  For  the  approximations  and  assumptions  made,  the  results  for  the  NIN  structures,  were  ordinary; 
the  departures  from  the  law  of  mass  action  were  negligible.  However,  for  the  NPN  structures,  the  results 
were  extraordinary.  Even  modest  values  of  P  (  near  10“/  m’)  resulted  in  two  orders  of  magnitude 
departures  from  the  mass  action  law.  These  results  which  are  quite  new  suggest  that  the  use  of  the  law  of 
mass  action,  one  of  the  most  frequently  used  relationships  in  the  design  of  devices  needs  to  be  re¬ 
examined.  This  issue  will  be  addressed  at  the  conference. 


This  work  is  supported  by  ONR  &  ARO. 


P9.  One-Dimensional  Analysis  of  Subthreshold  Characteristics  of  SOI-MOSFET  Considering 
Quantum  Mechanical  Effects,  Rimon  Ikeno,  Hiroshi  Ito  &  Kunihiro  Asada,  Dept,  of  Electronic  Eng., 
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Fac.  of  Engineering.  The  University  of  Tokyo,  7-3-1  Hongo,  Bunkyo-ku,  Tokyo  113,  Japan.  We  have  been 
studying  on  substrate  bias  dependence  of  subthreshold  characteristics  of  SOI  MOSFETs.  We  have 
simulated  subthreshold  characteristics  by  solving  simple  1-D  Poisson's  equation  on  SOI  multi-layer 
structure  (fig.l)  and  estimated  the  structural  parameters  of  real  devices.  Here,  we  considered  quantum 
mechanical  effects  in  electron  inversion  layer  of  thin  SOI  MOSIFET,  implementing  self-consistent  solver 
of  Poisson  and  Schrodinger  equation  in  the  1-D  subthreshold  simulator.  Figure  2  shows  an  example  of 
quantum  mechanical  results  of  band  diagrams  with  electron  distribution  probability  plots  in  SOI  layer  In 
Ais  case,  structural  device  parameters  are;  T,=7nm.T^,=30nm~-  7^  =  80nm  and  N,=1.5  x  10'’cm’. 
Figure  3  and  4  show  simulation  results  of  subthreshold  characteristics  of  SOI  MOSFET  for  2-D  1  D 
and  1-D  with  quantum  mechanical  effects.  Here  L  is  channel  length  and  is  substrate  bias  as  parameters. 
Classical  2-D  and  1-D  curves  show  good  agreement  in  subthreshold  region  in  long  devices,  though  they 
show  discrepancy  because  of  absence  of  drift  term  from  1-D  simulation  model.  Classical  and  quantum 
simulation  show  a  considerable  amount  of  difference,  which  proves  that  quantum  mechanical  analysis 
estimates  higher  threshold  voltage  than  the  classical  analysis.  We  calculated  V,,-V  characteristics  of  SOI 
MOSFET  (  V^(@Ij=10nA/L=W=  100pm),  with  classical  and  quantum  mechanical  simulator.  Figure  5 
shows  the  calculated  curves  with  measured  data.  Device  parameters  used  in  each  simulation  are  shown  in 
table  1,  with  designed  device  parameters  of  real  device.  Those  results  indicate  that  quantum  mechanical 
effects  must  be  considered  to  obtain  a  good  agreement  with  experimental  data,  based  on  the  reasonable 
device  parameters  in  thin  SOI  devices. 


10.  Total  Dielectric  Function  Approach  to  the  Electron  Boltzmann  Equation  for  Scattering  from  a 
Coupled  Mode  System,  B.A.  Sanborn,  Arizona  State  University,  Center  for  Solid  State  Electronics 
Research,  Tempe,  AZ  85287-6206.  The  total  dielectric  function  lends  itself  to  a  simple  and 

general  method  for  deriving  the  inelastic  collision  term  in  the  electron  Boltzmann  equation  for  scattering 
from  a  coupled  mode  system.[l]  Useful  applications  include  scattering  from  plasmon-polar  phonon  hybrid 
modes  in  modulabon  doped  semiconductor  structures.  In  the  Born  approximation,  the  inelastic  differential 
scattenng  rate  W”  can  be  expressed  in  terms  of  the  nonequilibrium)  er(g,0)),  including  both  electronic 
and  ionic  contributions.  Within  the  random-phase  approximation  W"'  separates  into  an  electron-electron 
interaction  and  an  electron-phonon  interaction  including  the  phonon  self-energy  due  to  polarization  of  the 
electrons.  'The  self-energy  modifies  phonon  dispersion,  so  that  hybrid  modes  appear  in  the  electron- 
phonon  collision  term.  An  iterative  method  will  be  deseribed  for  exactly  solving  the  Boltzmann  equation 
for  low  fields  including  dynamically  screened  electron  electron  scattering  and  plasmonphonon  coupling. 
Numencal  results  for  mobility  in  n-type  GaAs  show  that  inelastic  scattering  from  the  coupled  electron- 
phonon  system  dominates  over  acoustic  mode  scattering  at  77  K  for  moderate  doping,  if  mode  coupling 
and  electronelectron  scattering  are  included.[2]  New  results  to  be  presented  include  an  investigation  of 
electron  scattering  from  coupled  intrasubband-intersubband  electronic  excitations  in  modulation  doped 
quantum  wells.  ^ 


[1]  B.A.  Sanborn,  Phys.  Rev.  B  51, 14247  (1995). 

[2]  B.A.  Sanborn,  Phys.  Rev.  B  51,  14256  (1995). 


Pll.  A  Study  of  Transconductance  Degradation  in  HEMT  Using  a  Self-Consistent  Boltzmann- 
PoiMon-Schrodinger  Solver,  R.  Khoie,  Department  of  Electrical  and  Computer  Engineering,  University 
of  Nevada.  Ims  Vegas.  Las  Vegas.  NV  89154.  Recent  advances  in  III-V  compound  semiconductor  growth 
techniques  have  pushed  the  cutoff  frequency  of  high  electron  mobility  transistors  well  into  the  300  GHz 
ITt®  A  reported  design  and  fabrication  of  a  50-nm  self-aligned-gate  pseudomorphic 

AlInA^GalnAs  HEMT  with  a  maximum  transconductance  of  1740  mS/mm  and  a  cutoff  frequency  of  340 
GHz.  These  impressive  high-frequency  parameters  have  been  measured  at  very  low  gate  bias  voltages  In 
fact,  as  the  gate  bias  voltage  of  the  device  is  increased,  both  transconductance  and  the  cut-off  frequency  of 

^  have  been  reported  by 

others.[3][4]  The  high  transconductance  in  HEMT  devices  is  obtained  mainly  by  scaling  down  the  channel 
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length  of  the  device  to  500A'’.[5]-[6]  Further  reduction  in  the  gate  length  results  in  degradation  of  the 
transconductance  and  decrease  in  the  cut-off  frequency.  Kizilyalli,  et.  al.[7]  have  performed  a  Monte  Carlo 
study  of  short  channel  effects  in  a  submicrometer  AlGaAs/GaAs  HEMT,  and  suggested  that  the  high- 
frequency  performance  of  the  device  can  be  improved  by  scaling  the  gate  length  down  to  a  minimum  of 
about  0.1  |im,  beyond  which  the  device  transconductance  is  degraded  rather  rapidly.  They  discussed 
reasons  for  observed  degradation  of  the  transconductance  of  ultra-small  gate  length  HEMTs,  citing  poor 
charge  control  in  the  channel,  and  sharp  reduction  in  the  output  resistance  of  the  device  as  the  main  sources 
of  transconductance  degradation.  We  previously  reported  a  two-subband  self-consistent  Boltzmann- 
Poisson-Schrodinger  solver  for  high  electron  mobility  transistor  in  which  we  self-consistently  solved  the 
two  higher  moments  of  Boltzmann  equation,  along  with  Poisson  and  Schrodinger  equations.[8]  We  further 
incorporated  an  additional  self-consistency  by  calculating  field-dependent,  energy-dependent  intersubband 
and  intrasubband  scattering  rates  due  to  ionized  impurities  and  polar  optical  phonons.  In  this  work,  we 
have  used  our  Boltzmann-Poisson-Schrodinger  solver  and  studied  the  effects  of  the  intersubband  and 
intrasubband  scatterings  of  electrons,  on  the  transconductance  of  a  single  quantum  well  HEMT  device. 
The  results  of  our  simulations  (shown  in  Fig.  2)  exhibit  the  same  pattern  reported  by  [1]  and  [2].  The 
transconductance  increases  with  gate  voltage  at  low  gate  voltages,  and  then  decreases  as  the  gate  voltage  is 
further  increased.  Without  the  energy-dependent  field-dependent  scattering  rates,  the  transconductance 
decreases  linearly  with  increasing  gate  voltage.  We  concluded  that  the  degradation  of  transconductance  of 
the  device  with  applied  gate  bias  is  attributed  to  the  intersubband  and  intrasubband  scattering  mechanisms 
due  to  polar-optical  phonons  and  ionized  impurities. 


[1]  L.D.  Nguyen,  et  aL,  IEEE  Trans,  on  Electron  Devices,  VoL  39,  9,  pp.  1007-2014,  1992. 

[2]  Y.H.  Wang,  et  al.,  Solid  State  Electronics,  Vol.  37,  2,  pp.  237-241,  1994. 

[3]  R.  Plana,  et  al.,  IEEE  Trans,  on  Electron  Devices,  Vol.  40,  5,  pp.  852-858,  1993. 

[4]  S.K.  Sargood,  et  al.,  IEEE  J.  of  Quantum  Electronics,  Vol.  QE-29,  pp.  136-149, 1993. 

[5]  J.K.  Abrokwah,  et  al.,  IEEE  Trans,  on  Electron  Devices,  Vol.  40,  pp.  278-284, 1993. 

[6]  Y.  Zhao,  et  al.,  IEEE  Trans,  on  Electron  Devices,  Vol.  40,  1067-1070,  1993. 

[7]  I.C.  Kizilyalli,  et  al.,  IEEE  Trans,  on  Electron  Devices,  Vol.  40,  4,  pp.  234-249,  1993. 

[8]  R.  Khoie,  et  al..  Proceedings  of  the  International  Workshop  on  Computational  Electronics,  Edited  by  K. 
Hess,  U.  Ravioli,  and  R.  Dutton,  pp.  181-184,  1992. 


P12.  Advanced  Macroscopic  Device  Simulation  of  InGaAs/InP  Based  Avalanche  Photodiodes,  Joseph 
Parks  and  Kevin  F,  Brennan,  School  of  Electrical  and  Computer  Engineering,  Georgia  Tech,  Atlanta,  GA 
30332—0250;  and  Larry  Tarof  Bell  Northern  Research,  P,0,  Box  3511  Station  C,  Ottawa,  Canada, 
K1Y4H7.  In  this  paper,  we  analyze,  based  on  two-dimensional  drift-diffusion  and  hydrodynamic 
simulations,  how  variations  in  the  structural  components  of  an  InGaAs/InP  separate  absorption,  grading, 
charge,  and  multiplication  avalanche  photodiode  alter  the  device  performance.  The  model  is  employed  in 
conjunction  with  experimental  measurements  to  enhance  the  understanding  of  the  device  performance. 
Excellent  agreement  between  the  calculated  results  and  experimental  measurements  of  the  breakdown 
voltage,  dark  current,  mesa  punch-through  voltage,  photoresponse  and  gain  are  obtained.  We  examine 
specifically  how  variations  in  the  doping  concentrations,  doping  profiles,  and  layer  thicknesses  alter  the 
device  terminal  characterisitics.  In  this  way,  the  model  can  be  used  to  predict  how  processing  variations 
influence  the  overall  workings  of  advanced,  state-of-the-art  avalanche  photodiodes. 


PLl.  Self-Consistent  Scattering  Calculation  of  Resonant  Tunneling  Diode  Characteristics,  J,P,  Sun 
and  G.L  Haddad,  Center  for  High  Frequency  Microelectronics,  Dept  of  Electrical  Engineering  and 
Computer  Science,  University  of  Michigan,  Ann  Arbor,  ML  481 09-21 22.  We  perform  a  self-consistent 
calculation  of  resonant  tunneling  diode  (RTD)  I-V  characteristics  including  optical  phonon  scattering.  The 
self-consistency  is  botained  by  solving  the  Schrodinger  equation  and  Poisson’s  equation  iteratively  with  the 
Thomas-Fermi  approximation  used  for  the  device  contact  regions.  For  evaluation  of  phonon-assisted 
current  density,  the  optical  phonon  scattering  in  the  quantum  well  is  modeled  using  the  optical  model 
potential.  The  calculated  transmissions  and  electron  wavefunctions  illustrate  the  optical  model  and  effects 
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of  the  phonon  scattering  on  the  current  transport.  The  I-V  characteristics  we  obtain  from  the  model 
calculation  are  in  good  agreement  with  experimental  results.  This  work  manifests  the  importance  of 
including  both  self-consistency  and  optical  phonon  scattering  in  modeling  realistic  RTD  structures. 


PL2.  Microscopic  Theory  for  Transconductivity,  A.-P.  Jauho^M-  Bonsagef,  K.  Flensberg^^,  B,  Y,-K. 
//w",  and  J,  Kinaref,  Technical  University  of  Denmark,  DK-2800  Lyngby,  Denmark;  ^DFM,  A, 

Engelundsvej  I,  DK-2800  Lyngby,  Denmark  (permanent  address);  ^Chalmers  University  of  Technology, 
Goteborg,  S-41296  Sweden.  Measurements  of  momentum  transfer  between  two  closely  spaced  mesoscopic 
electronic  systems,  which  couple  via  Coulomb  interaction  but  where  tunneling  is  inhibited,  have  proven  to 
be  a  fruitful  method  of  extracting  information  about  interactions  in  mesoscopic  systems.[l]  Theoretical 
analyses  reported  so  far  have  been  based  on  the  semiclassical  Boltzmann  theory  [2],  memory  functions  [3], 
or  momentum  balance  equations.[4]  None  of  these  methods  possesses  the  full  rigor  of  a  diagrammatic 
Green  function  theory,  and  inclusion  of  higher  order  interaction  terms  necessarily  involves  some 
phenomenology.  We  report  a  fully  microscopic  theory  for  transconductivity  G12  ,  or,  equivalently, 
momentum  transfer  rate  between  the  system  constituents.  Our  main  formal  result  expresses  the 
transconductivity  in  terms  of  two  fluctuation  diagrams,  which  are  topologically  related,  but  not  equivalent 
to,  the  Azlamazov-Larkin  and  Maki-Thompson  diagrams  known  for  superconductivity.  We  have  evaluated 
the  general  expression  in  a  number  of  special  cases.  Our  new  results  include:  (i)  Semiclassical  clean  limit, 
when  Boltzmann  results  [2]  are  recovered;  (ii)  Dirty  systems  with  a  constant  relaxation  time,  which 
reproduces  the  memory  function  results  [3],  however,  for  energy  dependent  relaxation  times  the  final  result 
is  not  expressible  in  terms  of  standard  density-response  functions,  and  the  recently  reported  plasmon 
enhancement  of  the  drag  rate  [5]  reflects  the  details  of  the  underlying  scattering  mechanisms;  (iii)  At 
r  =  0  the  frequency  dependent  transfer  rate  can  be  computed;  (iv)  Weak  localization  corrections  can  be 
evaluated,  and  we  find  that  +Ag^^  ;  (v)  The  magnetic  field  dependence  of  (7^2 

computed,  and  we  can  prove  previous  conjectures  that  in  the  clean  limit  the  density-response  function  can 
be  replaced  by  its  magnetic-field  dependent  generalization.  Gi2iB)  is  strongly  enhanced  over  its  zero  field 
value,  and  it  displays  strong  features,  which  can  be  understood  in  terms  of  a  competition  between  density- 
of-states  and  screening  effects. 


[1]  T.-J.  Gramila,  et  al.,  Phys  Rev.  Lett.  Vol.  66,  1216  (1991). 

[2]  A.-P.  Jauho  and  H.  Smith,  Phys.  Rev.  B,  Vol.  47, 4420  (1993). 

[3]  L.  Zheng  and  A.H.  MacDonald,  Phys.  Rev.  B,  Vol.  48,  8203  (1993). 

[4]  H.C.  Tso,  et  al.,  Phys.  Rev.  Lett.  Vol.  68,  2516  (1992). 

[5]  K.  Flensberg  and  B.  Y.-K.  Hu,  Phys.  Rev.  Lett,  Vol.  73,  3572  (1994). 


P13.  Beating  in  the  RHEED  Intensity  Oscillations  During  Surfactant  Mediated  GaAs  Molecular 
Beam  Epitaxy:  Process  Physics  and  Modeling,  Vamsee  K.Pamula  and  R.Venkatasubramanian, 
Department  of  Electrical  and  Computer  Engineering,  University  of  Nevada,  Las  Vegas.  Las  Vegas,  NV 
89154.  In  a  recent  work,  it  was  shown  that  beating  in  the  reflection  high  energy  electron  diffraction 
(RHEED)  intensity  oscillations  can  be  obtained  during  MBE  growth  of  GaAs  when  the  surfactant,  Sn  is 
employed.[l]  It  is  also  observed  that  the  strength  of  beating  is  dependent  on  the  Sn  monolayer  coverage 
[1]  with  strong  beating  observed  for  0.4  monolayer  coverage  of  Sn.  Moreover,  the  period  of  oscillation 
decreases  with  Sn  coverage  even  though  the  actual  growth  rate  measured  by  high  resolution  transmission 
electron  spectroscopy  is  the  same.  The  origin  of  these  oscillations  and  its  dependence  on  Sn  coverage  have 
not  been  fully  explained.  In  this  work,  we  have  developed  a  rate  equation  model  of  growth  allowing 
elementary  kinetic  processes  such  as  adsorption  and  migration.  In  the  adsorption  model,  two  types  of 
possibilities  are  considered  depending  on  the  surface  of  adsorption,  i.e.,  GaAs  and  Sn.  In  addition  to  the 
normal  adsorption  process  considered  for  the  GaAs  surface,  the  Ga  atoms  arriving  on  the  Sn  surface  are 
allowed  to  incorporate  by  direct  exchange  process,  the  area  covered  by  the  Sn  is  assumed  to  grow  at  a 
slower  rate  compared  to  the  rest  of  the  area.  Results  of  RHEED  intensity  are  computed  as  a  function  of 
time  based  on  the  kinematic  theory  of  diffraction  assuming  that  the  electron  beams  reflected  from  the  Sn 
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covered  area  and  the  rest  of  the  area  are  incoherent.  This  is  justified  if  the  Sn  island  size  is  comparable  to 
the  coherent  length  of  the  electron  beam.  The  results  of  RHEED  intensity  oscillations  for  various 
monolaver  coverages  of  Sn  were  obtained  for  substrate  temperature  of  600°C  and  compared  with 
experimental  data.[l]  The  agreement  is  excellent  suggesting  that  the  proposed  model  of  growth  is  a 
plausible  one  for  Sn  mediated  GaAs  MBE  growth. 

[1]  G.S.Petrich,  A.M.Dabiran,  J.E.MacDonald  and  P.I.Cohen,  J.Vac.Sci.TechnoL,  B  9,  pp  2150-2153, 
(1991). 


P14.  Plasma  Process  Modelin  for  Integrated  Circuits  Manufacturing,  M.  Meyyappan  and  T.  R. 
Govindan,  Scientific  Research  Associates,  P.  0.  Box  1058,  Glastonbury,  CT  06033.  Plasma  processing  of 
semiconductors  is  widely  used  in  integrated  circuits  (IC)  manufacturing;  plasma  based  deposition  of 
dielectrics  and  etching  are  two  key  steps  in  device  fabrication.  At  the  microscopic  level,  prediction  of 
deposition  and  etch  profile  shapes  is  of  much  importance.  There  are  several  programs  such  as  EVOLVE 
and  SPEEDIE  to  address  these  needs.  These  programs  need  information  on  neutral  and  ionic  species  flux, 
ion  energy,  etc.  In  order  to  obtain  such  data,  one  needs  to  model  the  process  at  the  reactor  level.  The 
reactor  level  model  also  provides  other  valuable  information  on  deposition/etch  rates  and  uniformity  across 
the  wafer.  In  this  work  we  discuss  the  reactor  model  development  work  at  SRA.  The  reactor  model 
consists  of  discharge  physics,  chemistry,  fluid  flow  and  heat  transfer.  We  solve  the  governing  species 
conservation  equations  for  neutral  and  ionic  species,  flow  and  energy  equations  in  a  complex  manner. 
Homogenous  and  heterogenous  chemistry  are  included.  Our  code  includes  discharge  physics  and  power 
coupling  modules  suitable  for  rf  capacitively  coupled  reactors  and  high  density  ECR  and  ICP  reactors. 
Plasma  reactor  simulations  in  two-dimensions  are  computationally  intensive  due  to  the  coupling  of  physics 
and  chemistry  and  the  large  number  of  species.  Therefore,  we  have  also  developed  a  0-d  model  (well- 
mixed  reactor  model)  and  a  1-d  model  (in  the  flow  direction).  These  lower  order  models  take  much  less 
computer  lime  than  2-d  simulations  and  provide  valuable  information  (except  on  uniformity  issues)  for  use 
with  Eopographic  simulation  codes.  Results  from  these  three  codes  are  discussed  and  compared. 


MPl.  Derivation  of  an  Energy-Transport  Model  from  the  Boltzmann  Equation,  P.  Degond  and  N.  Ben 
Abdallah,  MIP,  (CNRS  UMR  9974),  University  of  Toulouse.  An  energy-transport  model  is  rigorously 
derived  from  the  Boltzmann  Transport  Equation  of  semiconductors  under  the  hypothesis  that  the  energy 
gain  or  loss  of  the  electrons  by  the  phonon  collisions  is  weak.  Retaining  at  leading  order  electron-electron 
collisions  and  elastic  collisions  (i.e.,  impurity  scattering  and  the  “elastic  part”  of  phonon  collisions),  a 
rigorous  diffusion  limit  of  the  Boltzmann  equation  can  be  carried  over,  which  leads  to  a  set  of  diffusion 
equations  for  the  electron  density  and  temperature.  The  derivation  is  given  in  both  the  degenerate  and  non¬ 
degenerate  case. 


MP2.  Topologically  Rectangular  Grids  in  the  ParaUel  Simulation  of  Semiconductor  Devices,  A. 

Asenov*,  A.R.  Brown,  S.  Roy,  and  J.R.  Barker,  Department  of  Electronics  and  Electrical  Engineering, 
University  of  Glasgow,  Glasgow,  G12  8QQ,  Scotland,  UK.  Topologically  rectangular  two  (2D)  and  three 
dimensional  (3D)  grids  are  particularly  attractive  for  the  parallel  simulation  of  semiconductor  devices, 
because  they  can  be  easily  partitioned  onto  processor  arrays,  enhancing  the  scalability  and  portability  of 
code.  Cartesian  finite  difference  grids  are  inherently  topologically  rectangular,  preserving  the  number  of 
nodes  along  the  horizontal  and  vertical  grid  lines.  However,  topologically  rectangular  finite  element  grids 
can  also  be  designed  for  a  wide  range  of  simulation  problems.  Fig.  1  represents  an  example  of  a  3D 
topologically  rectangular  grid  based  on  distorted  bricks,  used  for  the  triangulation  of  a  spherical  solution 
domain.  Such  a  grid  preserves  the  number  of  condensation  points  in  each  one  of  the  three  index  directions 
i  j,k  as  in  a  Cartesian  3D  grid.  The  spherical  solution  domain  can  be  mapped  into  a  cube  by  isoparametric 
transformations  mapping  each  of  the  distorted  bricks  into  unit  cubes.  There  are  significant  advantages  in 
speed-up  and  efficiency  when  such  topologically  rectangular  grids  are  partitioned  on  mesh  connected 
arrays  of  processors  instead  of  on  pipelines.  For  3D  grids  this  is  illustrated  in  Fig.  2  where  the  speed-up  of 
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a  linear  equation  solver  is  shown  as  a  function  of  problem  size  for  a  pipeline,  2D  and  3D  array  of 
processors.  This  figure,  based  on  detailed  performance  theory  [1],  also  shows  some  of  the  problems  arising 
from  the  'naive'  rectilinear  partitioning  of  the  grid.  Here  we  report  on  the  different  aspects  of  the 
implementation  of  topologically  rectangular  grids  in  the  simulation  of  semiconductor  devices.  As  an 
illustration  we  will  use  the  the  finite  element  3D  simulation  of  cellular  power  IGBTs.  The  schematic  view 
of  an  individual  cell  of  an  IGBT  is  illusU*ated  in  Fig.  3.  The  complex  shape  and  the  doping  distribution  in 
the  device  requires  a  3D  finite  element  discretization.  This,  together  with  the  computational  complexity  in 
the  simulation  of  power  devices,  makes  the  problem  a  distinguished  candidate  for  parallel  processing.  We 
use  spatial  decomposition  of  the  simulation  domain  over  arrays  of  processors.  The  basic  idea  of 
decomposing  a  3D  device  over  a  2D  array  of  processors  is  illustrated  in  Fig.4.  Compared  to  the  most 
frequently  used  pipeline,  the  use  of  a  2D  array  of  processors  minimises  the  inter-processor  communications 
by  reducing  the  ratio  between  the  bulk  and  the  surface  of  the  partition  subdomains.  A  universal 
communication  harness,  designed  in  PARDC,  provides  all  necessary  global  and  local  communications  and 
improves  the  portability  of  our  approach.  A  generalised  approach  to  the  generation  of  topologically 
rectangular  grids  has  been  developed.  In  the  IGBT  case,  because  of  the  symmetry  only  one  quarter  of  the 
device  has  been  discretized.  The  grid  conforms  to  the  shape  of  the  gate  and  the  metallurgical  pn  junctions. 
Fig.  5  illustrates  the  top  view  of  the  IGBT  grid  and  its  partition  over  a  3x3  array  of  processors.  A 
simulated  annealing  procedure  is  further  used  to  optimise  the  grid  partitioning.  The  final  partitioning  after 
the  annealing  is  illustrated  in  Fig.6.  For  parallel  matrix  generation  and  assembly  we  use  a  node  based 
procedure.  The  partition  sub-domain  is  scanned  node  by  node.  For  all  elements  which  condense  in  a  given 
node  only  the  contributions  to  this  particular  node  are  calculated.  This  leads  to  almost  100%  efficiency 
when  the  numbers  of  grid  points  are  equal  in  each  partition  sub-domain  (Fig.7).  A  Coloured  One  Step 
Block  Newton-SOR  scheme  has  been  adopted  for  the  Poisson  equation.  As  an  example  the  electric  field 
distribution  in  an  example  IGBT  simulation  is  given  in  Fig.  8.  The  applied  voltages  are  =  OV  and  = 
600V,  respectively. 


[1]  Asenov,  A.,  et  al.  Speed-up  of  scalable  iterative  linear  solvers  implemented  on  array  of  transputers. 
Parallel  Computing,  21,  (1995)  669-682. 


MP3.  Moment  Closure  Hierarchies  from  Kinetic  Theories,  C.  David  Levermore,  Department  of 
Mathematics,  University  of  Arizona,  Tucson,  AZ  8572 L  We  seek  models  that  properly  capture  the 
macroscopic  regime  when  the  mean  free  path  is  much  smaller  than  the  macroscopic  length  scales,  while  in 
the  transition  regime  give  values  for  physical  fluxes  (and  other  quantities)  that  are  hopefully  of  the  correct 
order  of  magnitude,  and  are  at  least  consistent  with  the  nonnegativity  of  the  particle  density.  By  doing  so, 
such  models  may  provide  a  bridge  over  the  transition  regime  that  may  be  useful  in  the  construction  of 
hybrid  macroscopic/kinetic  simulations.  This  paper  presents  a  systematic  nonperturbative  derivation  of  a 
hierarchy  of  closed  systems  of  moment  equations  corresponding  to  any  classical  kinetic  theory.  In  the 
context  of  gas  dynamics  the  first  member  of  the  hierarchy  is  the  Euler  system,  which  is  based  on 
Maxwellian  velocity  distributions,  while  the  second  member  is  based  on  non-isotropic  Gaussian  velocity 
distributions.  In  the  context  of  semiconductors,  one  obtains  generalizations  of  so-called  hydrodynamic 
models.  The  closure  proceeds  in  two  steps.  The  first  ensures  that  every  member  of  the  hierarchy  is 
hyperbolic  and  has  an  entropy.  This  implies  that,  unlike  for  many  traditional  approaches,  the  resulting 
equations  are  formally  well-posed.  The  second  step  involves  modifying  the  collisional  terms  so  that 
members  of  the  hierarchy  beyond  the  second  also  recover  the  correct  dissipative  behavior — the  Navier- 
Stokes  equations  for  gas  dynamics  and  the  drift  difusion  equation  for  semiconductors.  This  is  done 
through  the  introduction  of  generalized  BGK  type  collision  operators.  The  simplest  such  system  for  gas 
dynamics  in  three  spatial  dimensions  is  a  "M-moment"  closure,  which  also  recovers  the  behavior  of  the 
Grad  "1 3-moment"  system  when  the  velocity  distributions  lie  near  local  Maxwellians.  The  closure 
procedure  can  be  applied  to  a  general  class  of  kinetic  theories. 
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MP4.  Three-Dimensional  S-Matrix  Simulation  of  Single-Electron  Resonant  Tunnelling  through 
Random  Ionised  Donor  Statei^  Hiroshi  MizutajCentral  Research  Laboratory,  Hitachi,  Ltd,,  Kokubunji, 
Tokyo  185,  Japan,  Laterally-Confined  resonant  tunnelling  diodes  (LCRTDs)  exhibit  an  interesting 
interplay  phenomenon  of  3D  quantisation,  Coulomb  blockade  and  impurity-induced  tunnelling  depending 
on  the  structural  parameters.  This  work  specially  focuses  on  the  single-electron  resonant  tunnelling  (RT)  in 
LCRTDs  assisted  by  a  few  ionised  donors  in  an  unintentionally  doped  GaAs  quantum  dot.  The  3D  multi- 
mode  scattering  matrix  (S-matrix)  theory  [1]  is  adopted  for  this  purpose  newly  introducing  the  random 
scattering  potential  caused  by  ionised  donors  in  the  quantum  dot.  By  comparing  simulated  current-voltage 
characteristics  with  experimental  data  [2],  it  is  manifested  for  the  first  time  that  the  observed  current 
staircase  results  from  background  donors  and  so  reflects  the  random  configuration  of  the  donors  in  the 
quantum  dot.  In  the  present  simulation  the  3D  Schrodinger  equation  is  numerically  solved  with  scattering 
boundary  conditions  for  an  AlAs(4.2nm)/GaAs(5.6nm)/AIAs(4.2nm)  LCRTD  with  various  configurations 
of  a  few  ionised  donors.  The  S-matrix  and  hence  transmission  probability  is  calculated  based  on  the  2D 
lateral  eigenstates,  and  the  tunnelling  current  is  evaluated  in  the  global  coherent  tunnelling  scheme.  With  a 
few  ionised  donors  being  placed  in  the  dot,  the  calculated  energy-dependence  of  the  total  transmission  rate 
clearly  shows  new  peaks  and,  in  addition,  the  resonant  energies  are  dependent  on  the  donor  configuration. 
Visualised  electron  probability  density  reveals  that  these  resonances  originate  in  RT  via  single-donor- 
induced  localised  states.  This  new  type  of  RT  leads  to  current  steps  of  order  0.1  nA  near  the  current 
threshold  which  is  quantitatively  in  good  agreement  with  the  experimental  results. 


[1]  H.  Mizuta,  C.  Goodings,  M.  Wagner  and  S.  Ho,  J.  Phys.:  Condens.  Matter  4,  8783,  1992. 

[2]  C.  Goodings,  H.  Mizuta,  J.  R.  A.  Cleaver  and  H.  Ahmed,  J.  Appl.  Phys.  76,  1276,  1994. 


MP5.  Hierarchy  of  Full  Band  Structure  Models  for  Monte  Carlo  Simulation,  Umberto  Ravaioli, 
Beckman  Institute,  University  of  Illinois  at  UrbanaChampaign,  Urbana,  IL  61801,  Monte  Carlo  methods 
which  include  the  numerical  band  structure  of  semiconductor  materials,  provide  an  accurate  physical 
picture  of  hot  carrier  transport.  These  methods  find  direct  applications  in  transport  calculations  and  device 
simulation,  and  are  also  useful  for  interpretation  of  experiments  and  calibration  of  other  models  based  on 
simplified  band  structures.  Main  drawbacks  of  Monte  Carlo  approaches  are  the  computational  cost  and  the 
numerical  noise.  When  details  of  the  numerical  band  structure  are  utilized,  it  is  very  important  to 
understand  which  aspects  of  the  physics  are  relevant  in  order  to  design  efficient  simulation  approaches  that 
adequately  resolve  the  behavior  of  carrier  high  energy  tails.  This  talk  will  discuss  the  various  hiererachy 
levels  that  are  possible  when  the  full  band  structure  is  considered.  At  the  highest  level,  the  scatterings  are 
treated  using  complete  k-k'  transition  rates,  which  entail  extremely  memory  intensive  computational 
applications.  At  the  lowest  level,  the  scattering  anisotropy  is  neglected  and  the  scattering  rate  is  considered 
to  be  a  constant  average  value  on  energy  isosurfaces  of  the  bandstructure.  This  model  is  more  practical  for 
device  simulation,  and  examples  will  be  presented  which  show  how  reasonable  speed  and  efficiency  can  be 
obtained  for  practical  applications  on  workstations.  In  between  the  two  extremes,  it  is  possible  to  design 
intermediate  models  which  preserve  some  essential  features  of  both.  A  model  which  has  been  recently 
implemented  subdivides  the  momentum  space  into  a  set  of  equivalent  valleys,  using  as  a  guideline  the 
deformation  potential  variations  evaluated  inside  the  Brillouin  zone  from  k-k'  transition  calculations. 
Within  each  valley,  a  specific  deformation  potential  is  selected  for  the  phonon  scattering  rates,  and  the 
transport  is  treated  as  in  the  isotropic  case  for  more  efficient  computations.  At  all  levels  of  the  band 
structure  hierarchy  of  models,  there  are  similar  issues  of  numerical  noise,  related  to  the  sampling  of  real 
and  momentum  space  that  the  Monte  Carlo  method  necessarily  performs  with  a  relatively  small  number  of 
particles.  The  last  part  of  the  talk  will  address  computationally  efficient  approaches  based  on  the 
assignment  of  variable  weights  to  the  simulated  particles,  in  conjunction  with  careful  gather-scatter 
procedures  to  split  particles  of  large  weight  and  combine  particles  of  small  weight.  Such  methods  can  be 
designed  to  reduce  the  variance  of  statistical  noise  and  may  provide  better  efficiency  and  a  broader  range 
of  applications  when  full  band  structure  approaches  are  applied  to  probe  the  physics  of  hot  carrier  energy 
tails  and  rare  events  in  device  simulation. 


19 


P15.  Energy  Grid  Generation  for  Resolving  and  Integrating  Ultra-Fine  Resonances  in 
Quantum  Device  Simulation,  Gerhard  Klimeck\  Roger  Lake^  R.  Chris  Bowen\  Chenjing  L  Fernando\ 
and  William  R,  Frensley;  ""Eric  Jonsson  School  of  Engineering,  University  of  Texas  at  Dallas,  Richardson, 
TX  75083-0688;  ^Corporate  R&SD,  Texas  Instruments  Incorporated,  Dallas,  TX  75235.  We  describe  the 
crucial  algorithm  which  allows  us  to  efficiently  include  scattering  from  alloy  disorder,  interface  roughness, 
acoustic  phonons,  and  polar  optical  phonons  in  single  band  and  multi-band  quantum  device  simulations. 
Such  calculations  require  the  integration  of  electron  density  and  current  density  distributions  over  a  range 
of  1  eV  while,  at  the  same  time,  resolving  resonances  on  the  order  of  1  peV.  To  perform  this  task  using  a 
limited  number  of  energy  nodes  (*  200)  requires  an  extremely  well  optimized  energy  grid.  Typical  [1] 
adaptive  integration  schemes  place  nodes  by  trial  and  error,  checking  for  the  relative  error  due  to  the 
introduction  of  additional  nodes,  because  the  "location"  of  the  spectral  features  is  initially  unknown. 
Resonance  finding  algorithms  for  single  band  [2]  and  multi  band  [3]  tight  binding  models  have  been 
developed  and  implemented.  We  have  developed  a  general  energy  grid  generator  that  utilizes  the 
information  of  the  resonance  energies  and  life  times  of  the  states,  the  conduction  band  edges  and  the 
thermal  distribution  of  the  carriers  in  the  contacts.  The  grid  minimizes  the  absolute  error  in  the  integration 
of  the  spectral  features  avoiding  the  placement  of  nodes  in  energy  ranges  without  significant  integral 
contribution.  Also  the  grid  is  constructed  such  that  the  integral  contributions  from  segment  to  segment  are 
about  the  same  size  and  therefore  minimizes  the  numerical  round-off  errors.  We  will  describe  the 
algorithm  and  present  device  simulations  which  include  the  scattering  mechanisms  listed  above. 


[1]  W.  H.  Press,  S.  A.  Teukolsky,  W.  T.  Vetterling,  and  B.  P.  Flannery,  Numerical  Recipes  in  C,  The  Art  of 
Scientific  Computing,  2  ed.  (Cambridge  University  Press,  Cambridge,  1992). 

[2]  C.  L.  Fernando  and  W.  R.  Frensley,  J.  AppL  Phys.  76,  2881  (1994). 

[3]  R.  C.  Bowen,  W.  R.  Frensley,  G.  Klimeck,  and  R.  K.  Lake,  Phys.  Rev.  B  (1995). 


P16.  Gridding  and  Discretization  For  Divergence  Form  (Semiconductor-Like)  PDEs,  Donald  J.  Rose, 
Hai  Shao,  Ming-Yaq  Kao;  Department  of  Computer  Science,  Duke  University,  Durham,  NC  27708.  In  this 
presentation  we  will  discuss  gridding  and  discretization  for  a  PDE  (or  PDF  system)  of  the  form  (in 
2D)-V-j  +  /  =  0,  (1)  where  j  —  j  specified  as 

=r  =  1,2,  (2)  /  =  (3).  For  example,  j  could  take  the  specific 

form  j  =  aVw+  (4)  where  a  is  a  scalar  function  andjS  =  W  arises  in  the 

continuity  equations  of  the  drift-diffusion  semiconductor  equations;  here)3  =  ^E,  fl ,  a  mobility  and  E  the 

electric  field.  Our  approach  to  discretization  (in  2D)  will  be  to  use  the  “box  method"  (finite  volume  and 

finite  box  methods  using  Voronoi  boxes  at  a  set  of  appropriate  grid  points  in  the  region  specifying  the 

PDE.  By  using  the  divergence  theorem,  this  transforms  (1)  to  the  operational  equation 

.7 Jy  +  f  /  =  0,  (5)  whereB/  is  a  box  containing  grid  point  p,  with  boundary <?B..  After 
Jmi  J  B 

approximating  the  integrals  in  (5)  on  each  box  B- associated  with  the  unknown  u^,  we  obtain  the  coupled 

system  of  nonlinear  (or  linear)  algebraic  equations  representing  the  discretization.  The  transient  extension 

du 

of  equation  (1),  —  —  y*  7  +  /  =  0,  where  j  and /are  defined  as  in  (2)  and  (3),  is  easily  handled.  Hence 
dt 

most  of  our  exposition  in  the  presentation  will  deal  with  the  PDE  in  the  form  of  (l)-(3).  Our  intention  in 
the  presentation  will  be  to  consider  increasingly  more  general  cases  of  j  as  in  (2)  to  make  the  exposition 
tractable  and  oriented  to  an  engineering  viewpoint.  For  example,  the  case  of  (4)  for  a  and  being 
constants  is  interesting  and  instructive  even  in  ID.  The  case  when  a  in  (4)  is  replaced  by  a  2  x  2  (constant) 
matrix  motivates  much  of  the  discretization  algebra.  In  order  to  discretize  (5),  we  need  to  approximate  the 
integrals.  We  will  derive  a  discretization  based  on  the  constant  j  assumption,  namely,  that  j  is  constant 
pairwise  on  the  lines  of  the  Delaunay  triangles  dual  to  box  .  This  allows  us  to  approximate;  •  n  and  to 

assemble  the  parts  of  the  line  integral  which,  in  turn,  relates  the  variable  ,  to  the  other  corresponding 
to  Delaunay  neighbors,  say  Pi^ ,  of  the  box  mesh  point  p^.  We  will  show  that  the  constant;  discretization  is 


20 


a  natural  generalization  of  the  Scharfetter-Gummel  discretization  used  in  semiconductor  modeling.  When  j 
is  linear,  as  in  (4),  the  constant  j  discretization  produces  a  linear  system  to  solve  with  interesting  matrix 
properties.  We  will  discuss  these  matrix  properties  and  their  algorithmic  significance.  We  will  also 
examine  the  relation  between  such  matrix  properties  and  the  boxes  and  triangles.  The  generalization  to 
nonlinear  j  of  the  more  general  form  (2)  will  also  be  examined. 


P17.  A  New  Method  to  Recover  Vectorial  Electric  Fields  and  Current  Densities  from  Unstructured 
Meshes,  Daniel  C.  Kerr  and  Isaak  D.  Mayergoyz,  Dept  of  Electrical  Engineering,  University  of  Maryland, 
College  Park,  MD  20742,  In  the  numerical  simulation  of  semiconductor  devices,  the  vectorial  electric  field 
and  current  density  field  are  required  throughout  the  domain  of  simulation  to  compute  various  physical 
models,  such  as  mobility,  impact  ionization,  and  Joule  heating.  In  finite-box  simulations,  the  discretization 
and  solution  do  not  uniquely  define  the  fields  off  of  the  edges  joining  the  nodes.  The  fields  must  be 
reconstructed  from  the  projections  of  electric  field  or  current  density  along  the  edges.  The  existing 
recovery  methods,  the  method  of  comer  averages  [1],  and  the  method  of  least-squares  fitting,  are  ad  hoc 
means  to  average  the  data.  With  3-D  device  simulation  becoming  more  common,  a  mathematically  sound 
recovery  method  is  needed.  This  paper  describes  a  new  recovery  method,  which  is  based  on  the  edge 
elements  of  the  finite-element  theory.  The  method  of  edge-elements  directly  interpolates  vectorial  values 
defined  on  the  edges  of  an  element  into  the  interior  of  the  element.  The  vectorial  interpolant  f  of  the  edge 
values  Fij  into  the  interior  of  element  of  element  k  is  f  =  /  kF.e.. ,  where  e.  is  the  basis  function  of 

edge  ij  and  E*"  is  the  set  of  edges  of  the  element.  The  edge  basis  function  is  defined  in  terms  of  the  scalar 
finite-element  shape  function  and  d^  is  the  length  of  edge  ij.  Reconstruction  of  the  field  using  edge 
elements  yields  a  non-constant  vector  function  defined  on  the  element  which  has  the  following  properties: 
(1)  the  projection  of  the  edge-element  reconstruction  on  each  edge  reproduces  the  original  data,  and  (2)  the 
tangential  components  of  the  field  are  continuous  from  one  element  to  the  next.  A  simple  modification  of 
the  basis  functions  is  needed  when  the  edge-element  method  is  applied  to  non-simplex  elements.  In  this 
form,  edge  elements  are  suitable  for  2-  and  3-D  nonsimplex  elements.  The  original  and  modified  2-D 
edge-elements  in  standard  position  are  listed  in  Table  1.  This  element-averaged  J\  which  can  be  used  in 
device  simulators,  is  also  listed. 
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P18.  Parallel  Spectral  Multilevel  Solution  of  the  Augmented  Drift-Diffusion  Equations,  M,B,  Davis 
and  G,F,  Carey,  ASE  and  EM  Dept,  University  of  Texas  at  Austin,  Semiconductor  devices  in  the  sub¬ 
micron  range  exhibit  non-local  electric  field  effects  such  as  velocity  overshoot.  Standard  field-dependent 
mobility  models  in  the  drift-diffusion  equations  do  not  model  such  effects  since  calculation  of  the  mobility 
is  dependent  only  on  the  local  value  of  the  electric  field.  These  non-local  effects  can  be  captured  by 
augmenting  the  carrier  mobility  with  a  term  proportional  to  the  local  rate  of  change  of  the  electric  field. 
We  present  results  for  an  augmented  mobility  model  with  sensitivity  studies  for  the  augmentation 
coefficient.  Comparison  studies  to  other  augmented  mobility  models,  hydrodynamic  and  Monte  Carlo 
results  are  also  presented.  Spectral  (p)  finite  elements  are  used  to  discretize  the  equations  in  a  decoupled 
Gummel  algorithm,  and  the  multilevel  solution  strategy  uses  projections  between  bases  of  different  degree 
(level).  Hierarchic  bases  are  particularly  well  suited  since  the  element  matrices  and  vectors  are  nested  and 
the  projections  easily  defined  and  performed.  The  projection  methods  for  p-multilevel  are  particularly 
important  and  are  developed  and  analyzed  for  Lagrange  and  hierarchic  bases.  The  element-by-element 
(EBE)  parallelization  is  natural  for  the  finite  element  method,  and  if  basis  degree  is  used  to  specify  the 
multigrid  level,  an  EBE  strategy  is  natural  for  the  multilevel  technique  as  well.  Algorithm  scalability  and 
efficiency  is  analyzed  and  tested. 
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P19.  Distributed  Algorithms  for  Three-Dimensional  Semiconductor  Device  Simulations,  Mei-Kei 
leong  and  Ting~wei  Tang,  Department  of  Electrical  and  Computer  Engineering,  University'  of 
Massachusetts,  Amherst,  MA  01003.  The  parallel  simulation  of  semiconductor  devices  has  attracted  much 
attention  in  recent  years.  In  [1]  [2],  various  massively  parallel  algorithms  for  the  solution  of  linear 
equations  arising  from  three-dimensional  semiconductor  device  simulations  using  direct  and  iterative 
methods  have  been  studied.  In  [3],  stationary  iterative  methods  (such  as  Gauss-Seidel  and  Successive  Over 
Relaxation)  were  used  on  an  array  of  transputers  to  solve  two-dimensional  simulation  problems.  A 
speedup  of  14.3  on  16  transputers  has  been  reported.  Because  of  the  higher  performance-to-cost  ratio,  the 
coarse  grain  parallelism  is  a  very  good  alternative  for  semiconductor  device  simulations.  In  [4],  an 
optimum  speedup  of  three  on  a  network  of  four  workstations  for  a  two-dimensional  device  simulation 
using  various  domain  decomposition  methods  has  been  reported. 

In  this  work,  we  study  the  coarse  grained  parallelism  using  a  cluster  of  DECstations  at  the  University  of 
Massachusetts.  The  performance  of  two  iterative  methods,  JNGS  and  PBICGSTAB,  are  investigated.  The 
JNGS  method  for  approximating  the  solution  of  nonlinear  semiconductor  equations  is  naturally  parallel. 
The  nonlinear  problem  is  solved  in  a  point-wise  fashion  and  therefore  a  global  Jacobian  matrix  is  not 
necessary;  consequently  a  tremendous  saving  in  memory  can  be  expected.  However,  the  convergence  rate 
of  this  method  begins  to  deteriorate  as  the  number  of  domains  increases.  On  the  other  hand,  PBICGSTAB 
method  is  simply  a  parallel  version  of  the  BICGSTAB  algorithm  to  solve  the  linearlized  semiconductor 
equations.  This  method  together  with  a  good  preconditioner  is  robust  and  has  been  reported  to  be  capable 
for  solving  extremely  ill-conditioned  problems.  We  have  simulated  a  three-dimensional  submicron  N- 
MOSFET.  The  current  density  vector  shown  in  Fig.l  clearly  indicates  the  three  dimensional  feature  of  the 
result.  Table  1  shows  the  performance  of  the  JNGS  algorithm  applied  to  the  solution  of  the  3-D  MOSFET 
simulation,  with  gate  bias  of  2.0  volts  and  drain  bias  of  0.001  volts  on  a  cluster  of  DECstationSOOO.  The 
cluster  is  simultaneously  shared  by  other  users  and  some  fluctuation  in  the  performance  has  been  observed. 
Yet,  a  near  ideal  speedup  is  clearly  indicated.  Most  importantly,  very  large  problems  that  cannot  be 
handled  by  a  workstation  previously  can  now  be  solved  by  using  a  network  of  workstations  in  a  relatively 
short  time.  However,  as  the  drain  voltage  is  increased,  the  number  of  iterations  required  to  solve  the 
problem  also  increases.  The  deterioration  of  the  convergence  rate  is  due  to  a  larger  coupling  between 
Poisson’s  and  continuity  equations  at  high  drain  bias.  In  this  case,  the  BICGSTAB  method  is  preferred.  A 
comprehensive  comparison  between  the  two  methods  will  be  reported. 
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P20.  Quadrilateral  Finite  Element  Monte  Carlo  Simulation  of  Complex  Shape  Compound  FETs,  S. 
Babiker,  A.  Asenov\  J.  R.  Barker  and  S.  P.  Beaumont,  Nanoelectronics  Research  Centre,  Department  of 
Electronics  and  Electrical  Engineering,  University  of  Glasgow,  Glasgow  GI2-8QQ,  UK.  The  ensemble 
Monte  Carlo  simulation  approach  plays  a  particularly  important  role  in  the  simulation  of  compound  FETs 
where  both  steady-state  and  transient  device  behaviour  are  governed  by  velocity  overshoot  effects. 
However,  in  almost  all  ensemble  Monte  Carlo  studies  of  MESFETs  and  HEMTs,  planar  or  rectangular 
solution  domains  are  considered.  This  is  usually  dictated  by  the  rectangular  finite-difference  grid  used  for 
the  discretization  of  Poisson's  equation.  In  addition,  the  surface  potential  pinning  effects  are  normally 
neglected.  In  contrast,  all  modem  sub  micrometer  compound  FETs  are  single  or  double  recess  devices 
with  complicated  recess  geometry.  Very  often  the  length  of  the  free  recess  region  is  comparable  to  the 
length  of  the  gate.  Device  parasitics  like  access  resistances  and  coupling  capacitances  are  strongly 
dependent  on  the  shape  and  surface  conditions  of  the  recess.  These  parasitics  are  critical  to  the  device 
characteristics,  limiting  in  many  cases  the  device  performance  and  the  advantages  of  device 
miniaturization.  To  be  used  not  only  in  qualitative  research  but  in  the  practical  design  of  real  sub 
micrometer  FETs,  Monte  Carlo  simulation  codes  should  combine  the  extensive  transport  capabilities  with  a 
precise  description  of  the  device  geometry  and  proper  handling  of  surface  effects.  The  use  of  finite 
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elements  as  a  discretization  approach  offers  the  necessary  flexibility.  Here  we  report  on  a  new  Monte 
Carlo  (MC)  module  incorporated  in  our  Heterojunction  2D  Finite  element  FET  simulator  H2F  [1,2].  For 
the  first  time  this  module  combines  a  precise  finite-element  description  of  device  geometry,  with  realistic 
particle  simulation  of  the  non-equilibrium  hot  carrier  transport  in  ultra-short  recess  gale  compound  FETs. 
Quadrilateral  finite  elements  are  used  for  the  discretization  of  the  solution  domain.  The  grid  is  generated 
by  deformation  of  originally  rectangular  sub-domains.  The  Galerkin  method  with  linear  isoparametric 
mapping  has  been  adopted  to  approximate  Poisson's  equation.  An  efficient  bi-conjugale  gradient  method  is 
used  to  solve  the  resulting  system  of  equations.  The  motion  of  particles  within  the  solution  domain  is 
based  on  a  unified  modular  approach.  A  single  quadrilateral  element  is  the  building  block  of  the  MC 
module  where  the  trajectory  of  each  particle  is  traced,  both  in  momentum  and  real  spaces.  Each 
quadrilateral  finite  element  contains  a  single  material  and  has  a  uniform  doping  distribution  and  local 
scattering  table  attached  to  it.  The  boundaries  between  the  elements  may  be  homo-,  hetero-interface, 
contact  or  line  of  symmetry,  and  particles  reaching  the  boundary  are  transferred  or  reflected  according  to 
the  local  transition  probability.  Using  a  Single  Programme  Multiple  Data  (SPMD)  parallel  approach,  we 
are  running  H2F  and  its  MC  module  in  parallel  on  a  Parsytec  64  Supercluster  and  Parsytec  X'plorer.  This 
makes  it  possible  to  use  the  MC  simulations  for  practical  design  work,  generating  the  necessary  I-V 
characteristics  in  parallel.  The  capabilities  of  the  new  finite  element  MC  module  are  illustrated  in  example 
simulations  of  two  compound  FETs  fabricated  in  the  Nanoelectronics  Research  Centre  of  Glasgow 
University.  Fig.  1  represents  the  SEM  cross  sectional  view  and  the  H2F  simulation  domain  of  a  200  nm 
MESFET  with  55  nm  offset  between  the  gate  and  recess  edge.  Fig.  2  represents  the  SEM  cross  sectional 
view  and  the  H2F  simulation  domain  of  a  200  nm  gate  length  pseudomorphic  HEMT  fabricated  virtually 
without  offset  between  the  gate  and  the  recess  edges.  The  role  of  the  surface  potential  pinning  is  illustrated 
in  Fig.  2  a,b  where  the  longitudinal  component  of  the  electric  field  and  the  average  electron  velocity  in  the 
channel,  with  and  without  surface  potential  pinning,  are  plotted  for  the  MESFET  from  Fig.  1,  In  the 
transient  simulation  of  both  devices  terahertz  oscillations  were  observed.  The  switching  behaviour  of  the 
pseudomorphic  HEMT  is  illustrated  in  Fig.  4  for  different  doping  concentrations  in  the  cap.  The  Fourier 
transform  of  the  currents  from  Fig.  4  is  given  in  Fig.  5. 


[1]  A.  Asenov,  D.  Reid,  J.  Barker,  N.  Cameron,  S.  Beaumont,  Proc.  of  International  Workshop  on 
Computational  Electronics,  Leeds  University  Press,  45-49,  1993. 
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P21.  2D  Finite  Element  Method  Simulation  of  Lateral  Resonant  Tunneling  Devices,  Zhi-an  Shao, 
Wolfgang  Porod,  and  Craig  5.  Lent,  Department  of  Electrical  Engineering,  University  of  Notre  Dame, 
Notre  Dame,  IN  46556.  We  numerically  investigate  device  applications  of  lateral  resonant  tunneling 
structures  which  consist  of  a  transmission  channel  with  attached  resonators.  In  our  previous  studies,  we 
have  found  that  the  transmission  amplitude  possesses  resonance/anti-resonance  features  in  quantum 
waveguide  systems  with  resonantly-coupled  cavities.[l]  Here,  we  explore  the  utility  of  this  transmission 
feature  in  negative  differential  resistance  device  applications.  We  solve  the  2D  effective-mass  Schrodinger 
equation  with  current-carrying  boundary  condition  [2]  by  the  finite  element  method.  We  first  calculate  the 
transmission  probabilities  of  two  lateral  resonant  tunneling  devices  at  different  biases,  then  we  calculate  the 
1-V  characteristics  at  finite  temperature.  We  show  that  negative  differential  resistance  with  high  current 
peak-to-valley  ratio  can  be  engineered  in  2D  lateral  resonant  tunneling  structure  through  appropriate 
placement  of  the  transmission  resonance/anti-resonance  feature  on  the  real-energy  axis  (or  the  zero-pole 
pair  in  the  complex-energy  plane).  The  experimental  realization  of  these  two  devices  at  low  temperature  is 
also  discussed. 
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P22.A  Lattice  Boltzmann  Method  for  Extended  Hydrodynamical  Models  of  Electron  Transport  in 
Semiconductors^  S.Succi,  EXEC-IBM,  Roma,  Italy,  and  P.Vergari,  Dept,  of  Mathematics,  University  of 
Catania,  Italy.  Lattice  Gas  (LG)  and  Lactice  Boltzmann  (LB)  methods  have  been  much  in  use  in  the  recent 
years  to  simulate  a  wide  host  of  fluid-dynamic  flows,  ranging  from  creeping  flows  in  porous  media  and 
multiphase  flows  with  surface  tension,  to  homogeneous  incompressible  fully  developed  turbulent  flows,  to 
name  but  a  few.  Although  originally  targeted  to  'plain'  hydrodynamics  as  described  by  the  Navier-Stokes 
equations,  the  LB  method  lends  itself  to  a  number  of  remarkable  extensions  which  permit  us  to  handle 
generalized  hydrodynamic  flows  including  passive  scalar  transport,  Rayieigh-Benard  convection, 
magnetohydrodynamics  and  others.  In  this  paper,  we  shall  report  on  some  preliminary  efforts  to  construct 
an  LB  model  able  to  describe  the  generalized  hydrodynamic  equations  governing  the  dynamics  of 
relatively  hot  charge-carrier  in  semiconductor  devices.  The  basic  question  behind  this  effort  is  whether  a 
generalized  LB  scheme  can  be  devised  which  can  capture  the  higher  order  moments  of  the  ’true'  Boltzmann 
equation,  so  as  to  capture  physical  effects  beyond  the  reach  of  the  drift-diffusion  approximation. 
Preliminary  computations  indicate  that  indeed  such  a  model  can  be  found  under  appropriate  assumptions 
on  the  nature  of  the  local  thermal  equilibrium  of  the  charge-carriers.  Further  theoretical  and  numerical 
investigation  is  definitely  required  to  make  quantitative  assessments  of  the  present  model  versus  existing 
techniques. 


P23.  Second  Order  Newton  Iteration  Method  and  Its  Application  to  MOS  Compact  Modeling  and 
Circuit  Simulation,  Zhiping  Yu  and  Robert  W.  Dutton,  AEL  204,  Center  for  Integrated  Systems,  Stanford 
University,  Stanford,  CA  94305.  A  novel  second  order  Newton  iteration  scheme  has  been  developed. 
Different  from  other  higher  order  methods  [1]  it  not  only  exhibits  the  cubic  convergence  rate  but  also 
preserves  the  linear  (thus  robust)  nature  of  the  conventional  Newton  method.  The  second  order  derivative 
is  used  to  correct  the  slope  of  the  tangent  in  such  a  way  that  the  curvature  of  the  curve  is  naturally 
incorporated  in  the  iteration.  The  extra  computational  effort  in  evaluating  the  second  order  derivative  is 
partially  compensated  by  reduction  of  the  iteration  times,  and  more  importantly  by  providing  a  better  initial 
guess  when  the  function  parameters  are  changed.  The  applications  of  the  proposed  scheme  to  the  charge- 
sheet  MOS  compact  modeling  [2]  and  prediction  of  region  crossing  in  the  fast  circuit  simulation  [3]  are 
demonstrated.  It  is  expected  that  the  wide  recognition  of  this  technique  will  provide  a  vital  alternative  in 
solving  nonlinear  problems  and  in  particular  for  obtaining  an  approximate  solution  in  the  closed  form.  The 
algorithm  for  the  scheme  is  as  follows.  Given  a  system  of  nonlinear  equations,  f(x,p)  =  0, 

(1)  where  bold  symbols  represent  vectors  and  x  is  the  basic  variable(s)  such  as  the 
surface  potential  in  the  charge-sheet  model  and  p  is  the  parameter(s)  such  as  the  applied  bias.  Assuming 
for  now  A/  =  =  i,  the  conventional  Newton-Raphson  method  updates  the  solution  as  follows, 

X  •  —  X  Ajf',  Ax  =  -f{x  )/  {x)  (2)  where  i  is  the  iteration  count  and^  is  the  first  derivative. 


The  proposed  scheme  modifies  the  update  using  the  second  order  derivative,  in  the  following  way  (the 
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details  of  derivation  will  be  provided  in  a  full  length  paper)  Ax  =  —  — 

fx 


1  +  -  — / 


(3). 


It  essentially  changes  the  slope  of  the  tangent  to  better  reflect  the  curvature  of  the  function  at  solution  x . 
Because  Eq.  (3)  still  represents  a  straight  line,  the  scheme  is  robust.  Moreover,  this  scheme  can  readily  be 
applied  to  the  projection  of  the  initial  guess  when  p  changes.  The  following  expression  is  for  the  change  in 
jc  due  to  change  in  p  for  the  same  function  value. 


/,+4^  1^''"  (/,  +  4Ap)’ 
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In  view  of  the  critical  role  of  initial  guess  in  Newton  iteration,  the  extra  effort  is  worthwhile  to  ensure  the 
solution  be  found  with  the  iteration.  The  example  on  charge-sheet  MOS  modeling  shown  on  the  next  page 
confirms  the  merit  and  effectiveness  of  this  method.  The  iteration  times  can  be  cut  as  much  as  one  half. 


MP6.  Testing  Hydrodynamical  Models  on  the  Characteristics  of  One-Dimensional  Submicrometer 
Structure,  A  .Af.  Anile*,  0.  Muscato*,  S.Rinaudo#,  P.  Vergari*,  *Dip.  of  Mathematics,  University  of 
Catania,  Italy,  #S.G.5.  Thompson,  Catania,  Italy.  Recent  advances  in  technology  leads  to  increasing  high 
speed  performance  of  submicrometer  electron  devices  by  the  scaling  of  both  process  and  geometry.  In 
order  to  aid  the  design  of  these  devices  it  is  necessary  to  utilize  powerful  numerical  simulation  tools.  In  an 
industrial  environment  the  simulation  codes  based  on  the  Drift-Diffusion  models  have  been  widely  used. 
However  the  striking  dimension  of  the  devices  causes  the  Drift-Diffusion  based  simulators  to  become  less 
accurate.  Then  it  is  necessary  to  utilize  more  refined  models  (including  higher  order  moments  of  the 
distribution  function)  in  order  to  correctly  predict  the  behaviour  of  these  devices.  Short  of  a  direct  Monte 
Carlo  simulation  (which  requires  prohibitively  large  computational  cost  in  an  industrial  setting) 
hydrodynamical  models  have  been  considered  as  viable  simulation  tools.  Several  hydrodynamical  models 
have  been  considered  in  the  literature  with  various  degree  of  sophistication  and  completeness  ([2]-  [3]-  [5]- 

[6]-  [8]).  For  a  given  device  structure  the  various  hydrodynamical  models  can  give  widely  different  results 
for  the  velocity  and  energy  profiles  according  to  the  various  assumptions  made  in  the  model  (the  crucial 
parameter  seems  to  be  heat  conductivity  of  the  electron  gas).  However  the  different  velocity  or  energy 
profiles  are  hardly  accessible  to  experimental  detec  ion.  Therefore  it  is  mandatory  to  discriminate  among 
the  various  hydrodynamical  models  on  the  basis  of  their  results  on  the  output  characteristics  of  the  electron 
device  which  are  measurable  (I-V  curves).  We  have  analyzed  two  classes  of  hydrodynamical  models: 

•  bfields  hydrodynamical  models  and  hfields  drift-diffusion  model  [1]  a  general  purpose  two-dimensional 
device  simulator  developed  at  the  Bologna  University  to  perform  a  large  kind  of  analysis  like  steady-state, 
small  signal  and  transient.  It  resolves  the  semiconductor  equations  supporting  the  most  important  physical 
effects  using  a  large  range  of  models.  Default  parameters  are  used  both  in  the  drift-diffusion  version  and  in 
the  hydrodynamical  version  (Tp^ZO),  heat  conductivity  phenomenologically  represented  by  the 
Weidermann-Frau  law  with  c=-2.1  (arbitrarily)). 

•  Self-consistent  extended  hydrodynamical  models  with  relaxation  times  determined  from  Monte  Carlo 
simulations  (homogeneous  and  non-homogeneous).  Heat  conductivity  consistently  represented: 

-Linearized  Fourier  law  and  Non  Linear  (gradient  dependent)  heat  conductivity.  There  are  no  free 
parameters  because  all  of  them  are  consistently  determined  from  Monte  Carlo  simulations.  The  results  are 
shown  in  the  accompanying  figure.  We  plan  to  manufacture  an  electron  device  in  order  to  validate  the 
different  models. 
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MP7.  Monte  Carlo  Analysis  of  Anisotropy  in  the  Iransport  Relaxation  Times  for  the  Hydrodynamic 
Model,  R  Brunetti*  M.C.  VecchW,  and  M.  Rudan# ,  *Dipartimento  di  Fisica,  Universitd  di  Modena,  Via 
Campi  213/A,  1-41100  Modena,  Italy;  #Dipartimento  di  Elettronica,  Universitd  di  Bologna,  Viale 
Risorgimento  2,  1-40136,  Bologna,  Italy.  The  purpose  of  this  paper  is  an  investigation  of  the  anisotropy 
features  of  the  relaxation  times  used  in  the  hydrodynamic  model.  This  method  in  fact  considers  the 
moments  of  rank  0  through  3  of  the  Boltzmann  Transport  Equation  (BTE),  constituting  a  set  of  two 
second-order  nonlinear  PDEs  in  the  real  space.  Its  solution  provides,  in  addition  to  the  concentration  of 
carriers,  their  mean  velocity,  energy,  and  energy  flux.  In  the  derivation  of  the  model  it  is  possible  to 
incorporate  the  effects  of  the  collisions  into  a  set  of  generalized  relaxation  times  for  the  moments  of  rank  1, 
2,  and  3,  this  yielding  three  tensors  of  rank  2,  4,  and  6,  respectively.  The  rank  is  then  decreased  by  taking 
the  trace  of  the  corresponding  PDEs  [1],  thus  obtaining,  with  no  loss  of  generality,  the  following 
expressions  for  the  generalized  relaxation  times; 

=  T  J  n(w  -  w“‘),  (1) 

L  « V.  f  me(k)  =  T^'  nV .  (2) 

In  the  above, /is  the  distribution  function  and  n  =  f  fd'k  the  electron  concentration.  The  inverse  total 
scattering  rate  T  and  the  scattering-out  are  defined  by  l/r=J  H(k,k)d^k'  and 

~  ,k)d  k' .  Quantities  0),V,and^  are  the  mean  energy,  velocity,  and  energy  flux  of 

A-l  A-1 

the  electrons,  respectively.  In  practice,  tensors  Xp  and  Xp  are  supposed  to  reduce  to  scalars.  However,  in 
a  non-equilibrium  condition  the  existence  of  an  electric  field  breaking  the  cubic  symmetry  of  the  crystal  is 
expected  to  reflect  into  some  degree  of  anisotropy  of  the  distribution  function,  whose  moments  over  k 
would  then  acquire  a  tensor  nature.  It  is  therefore  of  interest  to  investigate  the  validity  of  the 
approximation  above,  and  to  this  purpose  it  is  necessary  to  determine  a  full  solution  of  the  BTE, 
incorporating  those  features  of  the  band  structure  that  are  relevant  to  determine  the  anisotropy  effects.  A 
Monte  Carlo  simulator  for  electron  transport  in  Si  [2]  has  been  used,  accounting  for  the  full  3D  electron 
dynamics  in  the  k  space  for  a  homogeneous  system  and  including  all  six  ellipsoidal  nonparabolic  valleys, 
as  described  in  [2].  The  generalized  relaxation  times  were  calculated  directly  from  Eqs. 

(1)  and  (2),  since  the  average  quantities  are  different  from  zero.  All  components  o£  and  X  have  also 

been  independently  evaluated  by  means  of  a  set  of  microscopic  correlation  functions,  as  described  in  [3]. 
This  procedure  allows  one  to  evaluate  also  the  tensor  components  involving  directions  orthogonal  to  the 
field  direction.  Numerical  results  were  obtained  at  7=  300  and  r=  77  A  as  functions  of  the  electric  field 
strength  and  direction. 
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MP8.  Monte  Carlo  Simulation  of  Hole  Transport  in  Strained  Si,^Ge.  and  of  Electron  Transport  in 
Strained  Si*, Gabriele  F.  Formicone  and  David  K.  Ferry,  Center  for  Solid  State  Electronics  Research. 
Arizona^  State  University,  Tentpe,  Az,  85287-6206.  .Monte  Carlo  simulations  are  used  to  study  the  transport 
properties  of  holes  in  a  sttained  Si,.,Ge,  layer  and  of  electrons  in  a  strained  Si  layer.  Mobility  enhancement 
ratio  (effective  mobility  in  strained  channel  over  effective  mobility  in  unstrained  channel)  is  extracted  for 
both  holes  and  electrons  and  compared  with  available  experimental  data.[l]  Using  an  electric  field  pattern 
computed  by  Medici  simulation  for  pMOS  and  nMOS  SiGe  based  devices  as  the  input  for  our  Monte  Carlo 
simulation,  we  also  study  the  carrier  drift  velocity  profile  vs.  the  distance  along  the  channel.  The  physical 
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model  employed  for  the  conduction  and  valence  bands  structure  in  strained  Si  and  strained  Si,  ,^Ge^  is  taken 
from  Rieger  and  Vogi.[2]  Surface  roughness  scattering  is  treated  according  to  Ferry. [3]  Preliminary 
results  show  that  electron  mobility  enhancement  ratio  saturates  at  a  Ge  concentration  equal  to  0.15-0.20. 
This  is  due  to  the  fact  that  the  energy  splitting  between  the  two  lowered  valleys  and  the  four  raised  valleys 
becomes  greater  than  the  intervalley  optical  phonon  energy,  suppressing  inlervalley  scattering  between 
them. 
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P24.  Time-Dependent  Solution  of  a  Full  Hydrodynamic  Model  Including  Convective  Terms  and 
Viscous  Effect,  Deyin  Xu  Ting-wei  Tang*  and  Sergei  S.  Kucherenko**;  *  Dept  of  Electrical  and  Computer 
Engineering,  University  of  Massachusetts,  Amherst,  MA  01003;  **  Dept  of  High  Energy  Density  Phvsics, 
Moscow  Engineering  Physics  Institute,  Moscow,  Russia.  Over  the  past  decade  the  hydrodynamic  transport 
model  has  been  extensivelv  used  in  the  simulation  of  submicron  semiconductor  devices  in  order  to  more 
accurately  predict  the  nonequilibrium  and  non-local  phenomena  occurring  in  these  devices.  Most  of  these 
simulations,  however,  are  limited  to  either  steady-state  solution  and/or  use  of  a  simplified  hydrodynamic 
transport  model.  Recently,  we  have  developed  a  svstem  of  hydrodynamic  transport  equations  based  on  the 
first  three  moments  of  the  Boltzmann  transport  equation  (BTE)  calibrated  by  Monte  Carlo  data[l].  This 
transport  model  differs  from  the  conventional  ones  in  that  the  first-order  moment  of  the  BTE  is  taken  w.r.t., 

velocity  (u)  rather  than  momentum  {tlk).  Thus,  the  average  velocity  V  is  treated  as  one  of  state 
variables  along  with  yA.n  and  IV  .  The  resulting  velocity  “transport”  equation  involves  a  second-order 

tenser,  A  =<  VV  >,  which  replaces  the  energy  tensor  U  —<  Vhk  >  in  the  conventional  momentum 

A 

transport  equation.  In  one  dimension,  A  is  modeled  as  [1]  A  =  W(1  + 

dV  W 

2.85(2aW)“77'^v - )  - ^  Cl(W) ,  where  the  third  term  represents  the  bulk  viscosity  and  T]  ==  2  [1],  [2]. 

dx 

To  solve  this  system  of  transport  equations,  we  used  the  second  upwind  for  spatial  discretization[3].  For 
time  stepping  scheme,  we  used  Bank's  TR-BDF2  composite  method[4].  As  an  application,  we  simulated 
various  1^ — structures,  one  of  which  is  shown  in  Fig.  1.  The  study  shows  that  the  steady-state  is 
established  within  a  few  picoseconds  after  all  external  voltage  is  suddenly  applied  at  r  =  0.  Figs.  2,  3,  4  and 
5  show  time  evolution  of  the  electric  field,  the  electron  velocity,  the  average  energy  and  the  current 
density,  respectively.  Note  that  the  electric  field  changes  suddenly  from  its  initial  value  to  a  constant  value 
of  E  =  33.3KV/cm  at  t  =  0+.  The  average  energy  increases  smoothly  with  time  from  its  thermal 
equilibrium  value  to  the  final  distribution  (Fig.  4).  However,  the  velocity  has  a  rapid  overshoot  (Fig.  3)  at 
around  0.1  ps  near  the  N"  —  N  junction  and  then  gradually  relaxes  to  a  much  lower  steady-state  value. 
About  the  same  time,  the  current  density  also  exhibits  a  similar  surge  (Fig.  5).  This  surge  of  current 
density  is  interpreted  as  onset  of  a  free  electron  plasma  oscillation.  After  approximately  one  quarter  of 
plasma  oscillation,  however,  the  relaxation  process  due  to  collisions  with  impurities,  photons,  etc. 
dominates.  In  Fig.  6,  the  steady-state  velocity  distribution  with  and  without  the  viscous  term  are  compared. 
As  expected,  the  viscous  term  has  a  dissipative  effect  and  substantially  suppresses  the  velocity  overshoot 
near  the  N*  —  N  junction  which  results  in  a  better  agreement  with  the  M.C.  data.  We  are  now  extending 
this  work  to  the  simulation  of  narrow-base  BJTs  in  which  the  viscous  effect  is  expected  to  play  a  more 
important  role. 
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P25.  Non-Parabolic  Inhomogeneous  Hydrodynamic  Models  for  Semiconductor  Device  Simulation, 

Arlvnn  W.  Smith  and  Kevin  F.  Brennan,  School  of  Electrical  and  Computer  Engineering,  Georgia  Institute 
of  Technology,  Atlanta,  GA  30332-0269.  In  this  paper,  we  discuss  different  hydrodynamic  model 
formulations  suitable  for  the  simulation  of  inhomogeneous  semiconductor  devices  without  using  the 
parabolic  band  approximation.  The  most  common  model  for  the  inclusion  of  non-parabolicity  utilizes  the 
Kane  dispersion  law  {tlkf'  12m  =  W(}  +  aW) .  An  alternative  choice  to  the  Kane  dispersion  model  is 
the  power  law  relation  of  Cassi  and  Ricco,  {tlkY  12m  =  xW^ .  Generalized  particle  and  energy  flux 
equations,  without  any  assumption  of  the  distribution  function,  are  presented  based  on  the  two  choices  of 
the  dispersion  relation.  The  physical  significance  of  the  band  non-parabolicity  is  discussed  as  well  as  the 
advantages/disadvantages  and  approximations  of  the  two  non-parabolic  models.  To  properly  account  for 
band  non-parabolicity  the  field  term  must  contain  a  non-parabolicity  factor  in  the  particle  flux  equation. 
This  is  true  for  both  non-parabolic  formulations  independent  of  the  distribution  function  chosen  and  has 
been  commonly  neglected.  Accounting  for  non-parabolicity  in  the  field  and  inhomogeneous  material  terms 
in  the  hydrodynamic  model  using  the  Kane  dispersion  law  severely  restricts  the  energy  range  over  which 
the  model  is  valid.  The  valid  energy  range  is  a  factor  of  2  to  3  times  smaller  than  is  commonly  assumed  for 
the  Kane  formulation  due  to  the  binomial  expansion  used  in  the  derivation  of  these  terms.  Finally,  closed 
form  particle  and  energy  flux  equations,  derived  assuming  a  heated  Fermi-Dirac  distribution,  are  presented 
in  a  form  suitable  for  implementation  in  device  simulators.  Extensions  of  the  model  to  more  accurately 
describe  both  the  distribution  function  and  band  structure  will  also  be  discussed. 


P26.  Hydrodynamic  Device  Modeling  with  Band  Nonparabolicity,  J.  Cai  and  H.L  Cu,  Department  of 
Physics  and  Engineering  Physics,  Stevens  Institute  of  Technology,  Hoboken,  New  Jersey  07030.  We 
present  a  semiconductor  device  model  based  on  a  set  of  quantum  mechanically  derived  hydrodynamic 
balance  equations  capable  of  dealing  with  arbitrary  electronic  energy  band  structures.  Compared  with  the 
hydrodynamic  balance  equations  for  parabolic  energy  bands,  the  main  difference  is  the  introduction  of  a 
variable,  ensemble-averaged,  inverse  effective  mass  tensor,  which  modifies  the  momentum  balance 
equation  directly  and  the  other  balance  equations  indirectly.  As  in  the  parabolic-band  case,  the  momentum 
and  energy  relaxation  rates  are  cast  in  the  form  of  electric  field  dependent  frictional  force  and  energy 
transfer  functions,  with  full  account  of  electron-electron  interaction  effects,  such  as  dynamical  screening, 
and  these  functions  are  calculated  using  the  electronic  dielectric  function,  which  is  normally  treated  within 
the  random  phase  approximation.  However,  these  calculations  are  now  modified  with  the  introduction  of 
the  ensemble-averaged  inverse  effective  mass  tensor  to  treat  the  arbitrary  energy  band  dependence  on  the 
crystal  momentum.  Device  modeling  issues  such  as  discretization  and  optimization  are  addressed  here. 
The  device  model  is  applied  to  a  semiconductor  with  a  Kane-type  conduction  band  to  demonstrate  the 
usefulness  of  the  new  approach.  We  compare  with  other  works  dealing  with  Kane  nonparabolicity  and 
whenever  possible,  show  explicitly  deviations  from  parabolic  expressions  and  results.  Numerical 
simulations  using  this  model  is  carried  out  for  a  submicron  ballistic  diode  structure,  using  the  Kane  model 
conduction  band,  and  the  results  are  compared  with  those  for  the  same  device  but  under  the  parabolic-band 
approximation. 


P27.  Computation  of  the  Spectral  Density  of  Noise  in  Bulk  Silicon  Based  on  the  Solution  of  the 
Boltzmann  Transport  Equation,  Alfredo  J.  Piazza  and  Can  E.  Korman,  Department  of  Electrical 
Engineering  and  Computer  Science,  The  George  Washington  University,  Washington,  DC  20052.  The 
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paper  presents  numerical  simulation  results  for  the  spectral  density  of  noise  due  to  current  fluctuations  in 
bulk  silicon  based  on  a  new  semiconductor  noise  model  [1],  The  mathematical  framework  of  the  noise 
model  is  the  machinery  of  stochastic  differential  equations  (SDE).  In  semiclassical  transport  theory,  the 
differential  equations  which  describe  the  motion  of  an  electron  in  a  semiconductor  can  be  interpreted  as 
stochastic  differential  equations  which  are  driven  by  inhomogeneous  randomly  weighted  Poisson 
processes.  These  processes  model  the  random  scattering  of  electrons  in  momentum  space  due  to  interband 
and  intraband  scattering.  The  solution  of  such  differential  equations  is  a  Markov  process  which  can  be 
characterized  by  a  transition  probability  density  function  and  satisfies  the  Kolmogorov-Feller  equation.  In 
the  case  of  semiclassical  transport,  this  equation  is  identical  to  the  linear  Boltzmann  transport  equation. 
Based  on  this  formalism,  the  key  computations  for  the  evaluation  of  the  autocovariance  function  of  current 
fluctuations  are  reduced  to  the  transient  solution  of  the  Boltzmann  transport  equation  (BTE)  with  special 
initial  conditions.  The  key  feature  which  differentiates  this  model  from  other  microscopic  noise  models  is 
that  this  approach  is  strictly  within  the  framework  of  semiclassical  transport.  Consequently,  this  approach 
directly  connects  the  noise  characteristics  with  the  properties  of  the  inhomogeneous  randomly  weighted 
Poisson  processes  which  describe  the  physics  of  scattering  in  the  semiclassical  transport  model.  The 
solution  method  for  the  solution  of  the  BTE  is  based  on  the  Legendre  polynomial  method.  The  stationary 
solution  is  represented  by  the  zero  and  first  order  Legendre  polynomials,  while  the  transient  solution  is 
represented  by  the  zero,  first  and  second  order  Legendre  polynomials,  due  to  the  particular  form  of  the 
initial  condition.  The  paper  describes  in  some  detail,  the  numerical  algorithm  for  the  solution  of  the 
transient  BTE  and  the  numerical  results  for  the  spectral  density  of  current  fluctuations.  The  numerical 
results  presented  in  this  paper  are  for  bulk  silicon  under  static  electric  fields  where  the  noise  is  due  to 
acoustic  and  optical  phonon  scattering. 


[1]  C.E.  Korman  and  I.D.  Mayergoyz,  "A  Semiconductor  Noise  Model  for  Semiclassical  Transport,"  in 
preparation. 


P28.  Hydrodynamic  Device  Modeling  with  New  State  Variables  Specially  Chosen  for  a  Block 
Gummel  Iterative  Approach,  Wenchao  Liang,  Daniel  C.  Kerr,  Neil  Goldsman,  and  Isaak  D.  Mayergoyz^ 
Dept,  of  Electrical  Engineering,  University  of  Maryland,  College  Park,  MD  20742,  We  report  on  a  new 
numerical  method  for  solving  the  HD  equations  which  is  specially  tailored  for  use  with  a  block-Gummel 
iterative  method.  The  motivation  for  using  the  Gummed  method  for  this  2D  simulator  is  to  establish  a 
basis  for  future  3D  calculations,  where  memory  limitations  inhibit  use  of  full  Newton  solvers.  With  this 
approach,  instead  of  using  standard  variables  { ^,n, },  we  define  new  state  variables 

{  These  new  variables  have  the  following  attributes:  (i)  The  new  variables  facillitate 


tailoring  the  HD  equations  specifically  for  use  with  a  Gummel  iterative  method;  (ii)  The  new  variables 
yield  transformed  HD  equations  so  that  the  left-hand  side  of  each  equation  is  linear  with  respect  to  a  new 
state  variable  with  which  this  equation  is  identified;  (iii)  The  transformed  equations  yield  well-conditioned 
coefficient  matrices  upon  discretization;  (iv)  The  new  variables  help  resolve  the  fast  spatial  variation  in 
carrier  concentration  and  carrier  temperature  which  occurs  in  semiconductor  devices.  These  attributes 
were  realized  by  developing  mappings,  from  the  standard  variables  to  the  new  ones,  which  facilitate 
expressing  each  of  the  HD  equations  in  forms  which  are  linear,  and  either  self-adjoint,  or  very  close  to  it. 
This  approach  yielded  the  following  relations  between  the  new  { (l>,u,Vyg„  ]  and  the  original  variables 


{ ^,np,Te  JP 


n  = 


^u\ 

*  K 


A  new  Scharfetter-Gummel-type  discretization  was  developed  for  each  of  the  transformed  linear  HD 
equations.  The  SG  discretization  yields  a  discrete  system  with  coefficient  matrices  which  are  well 
conditioned,  thereby  helping  to  avoid  numerical  problems  associated  with  their  solution.  The  coordinated 
use  of  the  new  variables,  the  Gummel  iterative  approach,  and  the  SG  discretization  yields  a  robust 
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approach  to  solving  the  HD  equations.  We  demonstrate  the  new  approach  by  applying  it  to  a  realistic  0.35 
pm  2-D  LDD  MOSFET  structure.  The  following  page  gives  some  example  results  including  the  predicted 
LDD  device  I-V  characteristics. 


P29.  Modelling  of  Hot  Acoustic  Phonon  Propagation  in  Two  Dimensional  Layers,  N,  A,  Bannov,  V.  V. 
Mitin,  and  F,  T,  Vasko,  Department  of  Electrical  and  Computer  Engineering,  Wayne  State  University, 
Detroit,  MI  48202.  Acoustic  modes  in  layered  structures  differ  substantially  from  modes  in  bulk  material. 
We  have  obtained  the  quantum  kinetic  equation  for  confined  acoustic  phonons  interacting  with  2D  electron 
gas  and  solved  it  for  the  case  of  phonon  transport  in  a  layered  structure.  We  use  the  Wigner  distribution 
function,  describe  the  phonon  syb-system;  here  m  is  a  discrete  index  for  mode  number,  q^ 

and  r,  are  the  in-plane  phonon  wave  vector  and  coordinate.  The  kinetic  equation  for  have  the 

fd .  . . 

to 


following  form 


a  d  \ 

—  II, r  ii)  zzT  (q  II, th),  where  S  ,  is  the  renormalised  due 

dvwj 


electron-phonon  interactions  phonon  group  velocity,  is  a  collision  integral  responsible  for 

decay  of  phonon  modes.  We  have  obtained  explicit  expressions  for  7^  (q  ii,r  ii)  and  for  the  case  of  a 


deformation  potential  interaction  of  2D  phonons  with  2D  electrons.  The  renormalization  of  the  acoustic 
wave  velocity  and  characteristic  decay  time  have  substantially  different  dependence  on  m  and  for  aq„«l 
and  for  aqj»l,  where  a  is  a  width  of  the  quantum  well.  We  have  numerically  solved  the  kinetic  equation 
for  the  thermal  pulse  propagation  along  a  layer.  Due  to  existence  of  several  confined  acoustic  modes  in 
layered  structures,  and  complicated  dispersion  relation  for  confined  phonons,  the  arrival  time  for  different 
phonons  varies  significantly.  This  fact  may  be  used  for  probing  electron-phonon  interaction  in  layered 
structures  by  time-of-flight  technique. 


P30.  Numerical  Simulation  of  Heat  Removal  from  Low  Dimensional  Nanostructures,  V.  V.  Mitin,  N.  A. 
Bannov,  R.  Mickevicius,  and  G.  Paulavicius ,  Department  of  Electrical  and  Computer  Engineering,  Wayne 
State  University,  Detroit,  MI  48202.  Heat  removal  from  active  elements  of  high  performance  ultralarge 
integrated  circuits  is  one  of  the  major  technical  problems  in  developing  new  generations  of  IC.  We  have 
investigated  the  acoustic  phonon  radiation  due  to  electron-acoustic  phonon  interaction  in  double  barrier 
quantum  wells  and  quantum  wires.  The  electron  and  phonon  density  matrices  are  governed  by  coupled 
kinetic  equations  which  have  been  solved  numerically  by  the  Monte  Carlo  technique.  The  effect  of  hot 
electrons  on  the  radiation  pattern  of  emitted  phonons  has  been  taken  into  account.  We  have  calculated  the 
angular  and  energy  spectrum  of  nonequilibrium  acoustic  phonons  radiated  from  low  dimensional 
structures.  A  major  peculiarity  of  acoustic  phonon  emission  is  a  violation  of  the  conservation  law  for  the 
phonon  wave  vector  component  which  is  perpendicular  to  the  direction  of  spatial  confinement.  Due  to  this 
violation  the  characteristic  energies  of  the  emitted  acoustic  phonons  greatly  exceed  the  energies  of  acoustic 
phonons  emitted  in  bulk  materials  and  are  order  of  several  meV.  Nevertheless,  such  phonons  propagate 
ballistically  over  macroscopic  distances  and  have  been  detected  experimentally.  Except  for  the  case  of  a 
strongly  degenerate  electron  gas,  the  along-structure  wave  vectors  of  acoustic  phonons  emitted  by  electrons 
in  nanostructures  are  smaller  than  the  transverse  wave  vectors.  Therefore  the  acoustic  phonons  carry 
energy  primary  in  the  direction  normal  to  the  nanostructure.  We  have  demonstrated  that  due  to  large 
acoustic-phonon  energies  and  strong  scattering  in  nanostructures  the  heat  removal  by  ballistic  fluxes  of 
acoustic  phonons  is  extremely  efficient. 


P31.  A  Hot-Hole  Transport  Model  Based  on  Spherical  Harmonics  Expansion  of  the  Anisotropic 
Bandstructure,  M.  Harrer,  H.  Kosina,  Institute  for  Microelectronics,  Vienna  Technical  University, 
Gusshausstrasse  27^29,  A-1040  Vienna,  Austria.  Monte  Carlo  transport  simulations  call  for  effective 
methods  to  calculate  the  free  flight  duration  and  to  choose  the  scattering  mechanism  and  the  state  after 
scattering.  We  propose  a  representation  of  the  valence  bands  using  an  expansion  into  a  series  of  spherical 
harmonics  that  is  capable  of  resolving  details  of  the  band  structure  both  at  the  center  and  at  the  boundary  of 
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the  Brillouin  zone.  The  basic  intention  of  our  method  is  to  simplify  the  calculation  of  the  integrated 
scattering  probability.  Assume  that  the  energy-wave-vector  relationship  is  given  in  polar  coordinates: 
£  =  E(k,Q).  Here,  £  is  the  energy  of  the  hole,  k  the  magnitude  of  the  wave  vector,  and  Q  denotes 
(0,0).  We  now  introduce  a  coordinate  transformation  (k,  Q)— >  (£,n)  ,  which  transforms  k  to  £  and 
lets  Q.  unchanged:  k  -  ¥i(£  For  any  given  the  function  K  is  defined  to  be  the  inverse  of  the 
function  £.  The  function  K  can  be  interpreted  to  describe  the  shape  of  an  equi-energy  surface  in  k-space. 
Inversion  of  a  function  is  possible  only  in  an  interval  where  the  function  is  monotonous.  By  inspection  of 
the  full  band  structure  one  finds  that  both  the  heavy  hole  and  the  split-off  bands  can  entirely  be  represented 
by  such  functions  K,  Above  a  hole  energy  of  Ex(3.04eV)  inversion  of  the  light  hole  band  is  no  longer 
unique.  In  this  work,  we  represent  the  function  A'  as  a  series  of  sperical  harmonics. 

K,  (£,a?  =  ■;^  X  S  (cose)r„(cos  0), 

/=0/w=0 

b^H.hSO,  (1). 

Derivation  of  the  scattering  rates  is  considerably  eased  by  taking  the  third  power  of  K  as  the  function  to  be 
expanded.  For  symmetry  reasons  non-vanishing  coefficients  only  exist  for  even  values  of  I  and  m  being  a 
multible  of  4.  With  (1)  a  set  of  functions  contains  the  whole  band  structure  information.  The 

essential  advantage  of  the  sperical  harmonic  expansion  of  the  valence  band  is  the  resulting  representations 
of  the  total  scattering  rates  and  of  the  distribution  of  the  scattering  angle.  Our  transport  model  accounts  for 
three  different  scattering  mechanisms,  namely  acoustic  deformation  potential  (ADP)  scattering  in  the 
elastic  approximation,  optical  deformation  potential  (ODP)  scattering  and  ionized  impurity  scattering 
(ION)  in  the  Brooks  and  Herring  formalism.  The  angular  distribution  functions  of  the  solid  angle  after 
scattering  are  also  given  by  spherical  harmonics  series.  The  rejection  technique  is  used  to  choose  the  after 
scattering  state.  The  free  flight  time  is  calculated  by  a  self-scattering  method.  The  functions  ^/„(£)  are 

represented  numerically  by  means  of  a  finite  elemente  method.  The  free  parameters  of  the  series  are 
determined  by  a  variational  approach.  The  number  of  harmonics  was  made  a  function  of  energy  ranging 
from  =  20  at  lower  energies  to  1^  =  60  at  higher  energies.  In  this  work  the  steady-state  hole  transport 
in  silicon  has  been  simulated  using  the  expansion  (1)  for  the  heavy  and  light  hole  bands  up  to  a  hole  energy 
of  =  3,04eV.  The  split-off  band  has  been  neglected.  Figure  (1)  shows  the  numerical  band  structure 
compared  with  the  series  representation.  The  numerical  band  structure  has  been  computed  by  a  nonlocal 
empirical  pseudopotential  method.  Figure  (2)  depicts  the  resulting  drift  velocities  in  comparison  to 
measured  data  and  Figure  (3)  the  simulated  average  hole  energy,  both  as  a  function  of  the  electric  field 
applied  in  characteristic  directions. 


P32.  An  Improved  Ionized  Impurity  Scattering  Model  for  Monte  Carlo  Calculations,  G.  Kaiblinger- 
Grujin  and  H.  Kosina,  Institute  for  Microelectronics,  Vienna  Technical  University,  Gusshausstrasse  27-29, 
A-I040  Vienna,  Austria,  We  have  developed  a  physically  based  ionized  impurity  scattering  model 
including  the  following  corrections  to  the  standard  Brooks-Herring  model.  First,  momentum  dependent 
screening  of  impuritites  by  conduction  electrons  is  taken  into  account  assuming  degenerate  statistics. 
Second,  the  effect  of  multi-ion-scattering  is  included.  Dynamical  screening  is  described  by  a  function  of 
both  the  tranferred  momentum  q  and  the  Fermi  level  [1].  Unfortunately,  this  function  is  represented  by  an 
integral  which  cannot  be  solved  analytically.  We  approximated  this  integral  by  an  analytical  expression 
which  has  exactly  the  same  behavior  as  the  original  integral  for  large  q  and  is  a  very  good  approximation 
for  small  q  for  arbitrary  degeneration.  The  advantage  of  this  approach  is  that  we  are  able  to  get  a  closed 
form  for  the  scattering  rate  without  changing  the  physics  of  the  underlying  problem.  With  higher  doping, 
the  average  distance  between  tango  impurities  becomes  smaller  and  the  neighboring  ion  potentials  overlap 
appreciably,  so  that  the  single-site-model  for  ionized  impurity  scattering  breaks  down.  Therefore  it  is 
necessary  to  consider  scattering  processes  at  two  ion  potentials  simultaneously.  Equally  charged  pairs  of 
impurities  scatter  up  to  twice  as  effectivly  than  monopoles  [2].  The  well-known  problem  of  very  large 
scattering  rates  at  small  angles  is  commonly  solved  by  using  a  method  after  Ridley  [3],  which  essentially 
cuts  off  the  scattering  rates  at  small  impact  parameters.  As  Ridley's  method  makes  the  agreement  of  theory 
and  experiment  even  worse,  we  generalized  this  method  in  that  we  allow  arbitrary  filter  functions.  The 
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filter  function,  which  cuts  off  large  values  of  the  impact  parameter,  is  chosen  such  that  the  mobility 
remains  unchanged.  Comparing  our  results  with  experimental  data  it  can  be  noticed  that  our  impurity 
scattering  model  improves  the  agreement  between  theory  and  experimental  data  significantly  and  that  it  is 
therefore  more  suitable  for  Monte  Carlo  calculations  than  the  classical  Brooks-Herring  model.  Despite  the 
greater  complexity  our  model  doesn't  consume  much  more  CPU  time.  The  new  model  represents  a  rather 
good  trade-off  between  an  exact  theory  and  an  accurate  approximation  that  is  applicable  for  simulation  of 
semiconductor  devices. 


[1]  W.-Y.  Chung  and  D.K.  Ferry,  “Dynamic  Screening  for  Ionized  Impurity  Scattering  in  Degenerate 
Semiconductors,”  Solid-State  Electron.,  Vol.  31,  no.  9,  pp.  1369-1374, 1988. 

[2]  J.  Meyer  and  F.  Bartoli,  “Effect  of  coherent  multi-ion  interference  on  ionized-impurity  scattering  in 
semiconductors,”  Physical  Review  B,  Vol.  30,  no.  2,  pp.  1026-1029,  1983. 

[3]  B.  Ridley,  “Reconciliation  of  the  Conwell-Weisskopf  and  Brooks-Herring  formulae  for  Charged- 
Impurity  Scattering  in  Semiconductors:  Third-Body  Interference,”  J.  Phys.  C:  Solid-State  Phys.,  Vol.  10, 
pp.  1589-1593,  1977. 


P33.  Simulation  of  Electron  Transport  in  Strained  Si/SiGe  Heterostructures,  Mahbub  Rashed.  W.-K, 
Shihy  S.  Jallepalli,  R*  Zaman,  T.J.T  Kwan*,  and  C,  Af.  Maziar,  Microelectronics  Research  Center,  The 
University  of  Texas  at  Austin;  *Los  Alamos  National  Laboratory,  New  Mexico,  Scaling  of  silicon 
MOSFET  channel  lengths  to  achieve  increases  in  speed  and  drive  current  has  pushed  the  limit  of  gate 
length  towards  0.1  pm.  An  alternate  approach  for  drive  current  enhancement  is  to  increase  the  mobility  of 
the  charge  carriers  in  the  channel.  Recent  reports  have  suggested  that  electron  mobility  is  enhanced  when 
electrons  flow  in  a  strained-Si  channel  pseudomorphically  grown  on  relaxed  (001)  Si,.^Ge^  [1,2,3].  The 
enhancement  of  mobility  in  strained-Si  is  due  to  both  the  suppression  of  intervalley  scattering  and  the 
lower  effective  mass  due  to  the  valley  splitting  of  the  six-fold  degeneracy  of  the  silicon  conduction  band 
minima.  This  lifting  of  degeneracy  results  in  an  upward  shift  of  the  four  transverse  valleys  and  lowering  of 
the  two  longitudinal  valleys  (in  energy).  In  this  work  electron  transport  in  strained  Si/Si,.^Ge, 
heterostructures  is  studied  using  a  both  a  bulk  Monte  Carlo  (MC)  tool  and  one  developed  for  2D  systems. 
The  MC  simulator  SLAPSHOT  [4]  has  been  modified  and  enhanced  to  investigate  electron  transport  in 
strained  Si/SiGe  systems.  The  bulk  MC  simulator  is  based  on  a  multiband  analytical  model  representing 
the  features  of  a  realistic  energy  bandstructure.  The  scattering  rate  computation  is  based  on  a  non-local 
pseudopotential  bandstructure.  Recent  results  of  semi-empirical  pseudopotential  bandstructure  calculations 
show  that,  to  first  order,  biaxial  strain  yields  changes  in  only  the  relative  energy  and  not  the  shape  of  the 
valleys  [3].  Therefore,  in  our  work,  effective  masses  and  nonparabolicities  of  different  valleys  are  assumed 
to  be  unaffected  by  the  strain.  Fig.  1  shows  the  velocity-field  characteristics  of  strained  and  unstrained 
silicon  for  different  directions  of  the  field.  Figs.  2  and  3  illustrate  the  velocity-field  characteristics  at  300  K 
and  77  K  for  several  mole  fractions  of  relaxed  Sij.^Ge^.  Low  field  mobility  is  enhanced  by  about  80%  at 
300  K  and  by  35%  at  77  K,  as  compared  to  bulk  silicon.  The  higher  mobility  observed  at  room 
temperature  is  due  to  greater  suppression  of  intervalley  phonons  as  phonon  scattering  increases  with 
temperature.  While  bulk  simulation  ignores  the  quantization  effect  in  the  transport,  the  nature  of  transport 
in  the  inversion  layer  requires  consideration  of  these  effects.  The  subband  structure  is  calculated  by 
applying  the  formalism  based  on  the  effectivemass  approximation  with  bulk  non-parabolic  E(K)  relation, 
as  described  in  detail  in  [5]  for  the  case  of  Si  MOSFETs.  Scattering  of  the  2D  electrons  due  to  bulk 
phonons  and  surface  roughness  is  included.  The  formalism  developed  by  Price  [6]  has  been  adopted  for  the 
phonon  scattering  rate  calculation.  Ionized  impurity  scattering,  being  important  only  in  the  presence  of 
large  interface  charge  density  or  when  the  channel  is  weakly  inverted,  is  ignored.  Surface  roughness 
scattering  is  simply  modeled  with  the  formalism  of  Cheng  et.  al.[7],  A  long  channel  strained  Si  nMOS 
structure  with  uniform  substrate  doping  concentration  and  zero  source-drain  bias  is  simulated.  At  each 
gate  bias,  a  two  dimensional  potential  profile  for  the  device  is  obtained  by  solving  drift-diffusion  (DD)  and 
Poisson  equations.  Since  a  uniform  inversion  layer  is  assumed,  a  slice  along  the  depth  is  selected  at  the 
middle  of  the  channel.  One  dimensional  Schrodinger  and  Poisson  equations  are  then  solved  self- 
consistantly  for  the  slice  in  order  to  calculate  subband  structures  in  the  strained  silicon  channel.  This 
calculation  is  followed  by  a  single  particle  MC  simulation,  where  a  single  particle  is  launched  according  to 
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the  density-of-states  (DOS)  weighted  equilibrium  distribution.  Fig.  4  shows  the  electron  concentrations  in 
the  strained  channel  for  different  gate  bias  conditions.  The  enhancement  of  effective  electron  mobility  in 
strained  silicon  channel  nMOS  structures,  as  compared  to  conventional  silicon  nMOSFETs,  is  presented. 
Fig.  5  shows  the  improved  mobility  in  strained  Si  nMOS  as  compared  to  silicon  nMOS  structure. 
Agreement  between  calculated  and  experimental  [9]  enhancement  of  effective  mobility  is  shown  in  Fig  6. 
This  agreement  is  sufficiently  encouraging  that  this  tool  will  be  used  for  preliminarv  evaluation  and 
analvsis  of  Si/Si, .,Ge,  based  devices. 
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P34.  A  New  Concept  for  Solving  The  Boltzmann  Transport  Equation  in  Ultra-fast  Transient 
Situations,  Ming-C.  Cheng,  Department  of  Electrical  Engineering,  University  of  New  Orleans,  New 
Orleans,  Louisiana  70248.  It  was  proposed  that  evolution  scales  of  f(k),  which  obeys  the  Boltzmann 
transport  equation,  can  be  characterized  by  relaxation  times.  This  concept  illustrated  in  Fig.  1  is  based  on 
the  fact  that  the  kinetic  distribution  function,  f(k),  can  be  described  by  the  infinite  set  of  moments  (i.e., 
hydrodynamic  parameters);  or/(/:)=  f (k^n^k^e^K  The  axis  in  Fig.  1  represents  the 

temporal  scale  of  evolution  for  corresponding  to  the  scales  of  moments.  In  semiconductor, 
>  Tg  >  ,  and  the  characteristic  times  of  higher-order  moments  are  assumed  to  be  smaller  than  . 

Information  infik),  after  a  drastic  change  in  field,  described  by  the  smaller-scale  moments  tends  to  vanish 
in  a  shorter  time.  As  shown  in  Fig.  l,f(k)  tends  to  evolve  through  f^(kyn,£,k)  (a  -scale  hydro- 
kinetic  distribution)  and  (a  -scale  hydro-kinetic  distribution),  and  eventually  into  a  quasi¬ 

equilibrium  distribution,  fgik,n).  In  fast  transient  situations,  fe  or  f^  can  be  chosen  to  describe/  The 
approach  to  fs  has  been  introduced  in  Ref.  1.  In  the  present  study,  evolution  of  the  hydro-kinetic 
distribution  fromf^  to  fc  is  assumed  to  be  a  relaxation  process  described  by  S^k,  and  the  change  in  field. 
The  process  can  be  performed  numerically.  Results  for  electrons  in  Si  <  100>  subjected  to  an  ultra-fast 
step  field  are  illustrated  in  Figs.2-4  where  a  two-valley  model  is  used.  Figs.  3  and  4  clearly  show  that  fs 
responds  to  the  change  in  field  more  slowly  than /and/,  during  the  velocity  overshoot  interval  where  the 
distribution  function  is  strongly  influenced  by  the  velocity  relaxation.  The  velocity  dependence  is  not 
included  in  the  -scale  distribution,  /g.  Inclusion  of  the  velocity  relaxation  in  /  significantly  improves 

the  accuracy  of  the  hydro-kinetic  distribution,  as  illustrated  in  Figs.  3  and  4  where/ and  /  evolve  closely. 
Calculations  for  the  evolution  take  only  about  10  seconds  on  a  486/33  PC  for  a  single-valley  band  model 
and  25  seconds  for  a  2- valley  model.  The  computational  model  for  the  evolution  will  be  presented. 

[1]  Ming-C.  Cheng  and  Rambabu  Chennupati,  J.  Phys.  D,  28,  160  (1995). 


TAl.  Recent  Advances  in  Device  Simulation  Using  Standard  Transport  Models,  Giorgio  Baccarani, 
Massimo  Rudan  and  Martino  Lorenz\ni,  Dipartimento  di  Elettronica,  Informatica  e  Sistemistica,  Viale 
Risorgimento,  2-40136  Bologna,  Italy.  Numerical  simulation  of  semiconductor  devices  has  reached  a 
mature  development  stage.  Two-and  three-dimensional  simulation  codes  are  increasingly  being  used  both 
for  a  better  understanding  of  device  behavior  and  for  design  purposes.  Furthermore,  commercial  tools 
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featuring  a  friendly  user  interface  and  sophisticated  simulation  environments  are  currently  made  available 
to  interested  users.  Yet,  it  is  fair  to  say  that  a  number  of  unsolved  practical  and  theoretical  problems  still 
exist  which  make  it  difficult  to  achieve  an  accurate  prediction  of  the  device  performance.  Such  problems 
range  from  structure  definition  to  mesh  generation  in  3-D;  also,  discretization  techniques  and  numerical 
stability  are  still  far  from  satisfactory  and  some  specific  physical  models  turn  out  to  be  very  hard  to 
quantitatively  predict.  In  this  paper,  we  address  some  the  above  problems  from  an  engineering  standpoint 
and,  whenever  possible,  indicate  desirable  research  areas  waiting  for  new  solutions.  Some  of  the  above 
limitations  are  listed  below. 

•  When  devices  are  inherently  three-dimensional,  a  realistic  definition  of  the  device  structure  turns  out  to 
be  very  hard  to  achieve  with  the  necessary  accuracy  in  both  device  morphology  and  impurity  profiles; 

•  Automatic  mesh  generation  in  3-D  is  still  a  very  challenging  task.  This  is  due  to  the  constraints  imposed 
by  the  sensitivity  of  the  numerical  solutions  to  the  quality  and  size  of  the  geometrical  elements,  and  by  the 
need  of  heavy  mesh  refinement  whenever  rapidly-varying  functions  are  to  be  accurately  determined; 

•  The  most  popular  discretization  scheme,  which  is  based  on  finite  boxes,  is  based  on  the  assumption  of  a 
spatially  uniform  current  density  along  the  element  sides.  Such  an  assumption  is  actually  correct  only  when 
the  current  density  is  parallel  (or  perpendicular)  to  the  element  side,  but  fails  when  the  current  flows  in  a 
slanting  direction  with  respect  to  them.  The  problem  is  especially  severe  at  the  drain  end  of  the  channel  in 
MOSFETs,  where  impact  ionization  occurs  under  saturation  conditions; 

•  The  electric  field  is  known  with  poorer  accuracy  than  the  electric  potential,  as  the  former  is  basically 
obtained  as  a  difference  between  close  potential  values.  As  a  result,  all  field-dependent  physical  effects, 
such  as  impact  ionization,  band-to-band  and  trap-assisted  tunneling,  Pool-Frenkel  enhanced  emission  rates 
and  Fowler-Nordheim  tunneling  across  the  oxide  turn  out  to  be  affected  by  relatively-large  errors; 

•  Hot  carrier  injection  into  the  gate  oxide,  which  is  the  dominant  mechanism  for  the  programming  of 
EPROM  and  flash  EEPROM  cells,  is  very  hard  to  reliably  predict  as  the  technology  changes.  This  is  due 
to  the  extreme  sensitivity  of  the  energy  distribution  function  at  large  energies  upon  the  physical  model 
(band  structure  and  scattering  mechanisms). 

This  paper  addresses  the  above  limitations,  and  discusses  new  solutions  and  perspectives  to  overcome 
them. 


TA2.  An  Efficient  Solution  Scheme  for  the  Spherical  Harmonics  Expansion  of  the  Boltzmann 
Transport  Equation  Applied  to  Two-Dimensional  Devices,  Maria  Christina  Vecchi\  Jan  Mohring  ^  and 
Massimo  Rudan\  ^  Dipartimento  di  Elettronica,  Universitd  di  Bologna,  Viale  Risorgimento  2,  40136 
Bologna,  Italy.  ^  Universitdt  Kaiserslautern,  FB  Mathematik,  E.  Schrodingerstrasse,  D’67663 
Kaiserslauter,  Germany.  In  recent  years  the  SHE  of  the  BTE  has  been  successfully  tested  in  a  wide  range 
of  problems  in  the  field  of  electron  transport  simulation  [1,  2].  Thanks  to  the  assumption  of  spherically- 
symmetric  bands,  the  differential  equation  in  energv  and  space  can  be  transformed  into  a  differential 
equation  in  space,  and  a  difference  equation  in  energy  [3].  This  allows  one  to  connect  to  each  point  in 
space  only  those  nodes  in  energy  coupled  via  the  scattering  operator.  On  the  other  hand,  a  detailed 
knowledge  of  the  distribution  function  at  low  energies  is  necessary  to  correctly  compute  its  averages 
(namely  concentration,  mean  energy,  mean  velocity,  etc.).  Because  of  this,  several  solutions  on  grids 
displaced  in  the  energy  domain  have  to  be  carried  out.  We  propose  here  a  new  scheme  that  improves  the 
solution  at  low  energies,  keeping  the  desired  accuracy  in  the  calculation  of  the  mean  quantities  while 
saving  a  significant  amount  of  CPU-time.  This  is  important  in  view  of  the  applications  of  the  method,  since 
the  typical  number  of  nodes  to  be  used  in  the  combined  space-energy  domain  is  in  the  range  of  10^-10^ 
Solution  Method:  The  SHE  method  in  steady  state  provides  the  second-order  difference-differential 
equation  [3] 

+  3c^(s*[w;/;  -  -  who'll  =  0,  (1) 

in  the  unknown  distribution  function  /q  .  The  equation  is  solved  in  a  three-dimensional  domain  (x,  y,  H) 
with  boundaries  defined  in  space  by  the  boundaries  of  the  device  and  in  energy  by  = — Q(pix,y) )  and 

— <l(p(x,y) ,  where  is  the  maximum  energy  of  the  band  system.  A  prismatic  mesh  with 


dr, 


dr, 
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triangular  elements  in  the  (jc,  y)  plane  is  adopted.  The  nodes  at  constant  total  energy  H  are  uniformly 
spaced  in  energy  by  intervals  A//  =  In,  with  n  integer.  Thanks  to  the  absence  in  (1)  of  partial 

derivatives  with  respect  to  H  each  node  is  connected  along  the  H  direction  only  to  the  nodes  at 
(x,  J,  //  ±  )  via  the  phonon  scattering  operator.  The  resulting  algebraic  system  is  thus  decomposed 

in  n  decoupled  subsystems.  In  the  solution  scheme  we  present  here,  eqn.  (1)  is  solved  in  the  full  energy 
domain  H^]  only  once,  while  several  additional  solutions  are  computed  in  the  reduced  energy 
domain  [H^^,  fnxHJ,  where  =flO)„p  ^q(pix,y)»  The  terms  in  eqn.  (1  )  which  are  external  to  m 

X  are  expressed  by  interpolating  the  values  previously  computed  on  the  full  domain.  By  letting  m  = 
10  one  can  express  the  values  of  the  distribution  function  outside  the  domain  by  an  interpolation  scheme 
that  maintains  the  linearity  of  the  system;  conversely,  to  compute  the  local  solution  with  m  =  1  one  has  to 
resort  to  a  more  sophisticated  interpolation  scheme.  It  should  be  added  that,  in  addition  to  the  loss  of 
linearity  of  the  system,  a  local  solution  scheme  with  m  =  1  does  not  allow  one  to  correctly  compute  the 
currents  in  the  device  because  of  the  early  truncation  of  the  solution  domain.  We  present  here,  for  the  case 
m  =  10,  the  comparison  of  the  results  of  this  solution  scheme  (consisting  of  1  full  solution  and  9  local 
solutions)  with  those  of  a  full  solution  computed  on  10  grids  in  The  test  device  is  a  2D-MOS 

transistor  with  effective  channel  length  0.25mm.  Figs.  1  and  2  present  the  comparison  of  the  electron 
concentrations  and  the  normalized  electron  mean  energies  along  a  section  of  the  device  parallel  to  the  si- 
siOj  interface,  and  Fig.  3  compares  the  drain  currents.  For  comparison  the  current  computed  by  HYDRO- 
HFIELDS  [4],  using  a  mobility  model  consistent  with  SHE,  is  also  reported  in  Fig.  3.  In  conclusion,  the 
scheme  based  on  the  solution  of  n  subsystems  displaced  in  the  full-energy  domain  can  be  replaced  without 
loss  of  accuracy  by  a  scheme  which  considers  a  single  complete  solution  and  (n — 1)  local  solutions  in  the 
interval  m  x  Thanks  to  its  structure,  the  new  method  is  intrinsically  more  efficient.  In  the 
example  shown  here,  a  reduction  factor  of  about  5  in  the  CPU-time  has  been  achieved. 
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TA3.  Brownian  Approach  to  the  Simulation  of  Small  Semiconductor  Device^  C.R.  Arokianathan,  A, 
Asenov"^,  and  J,H.  Davies,  Nanoelectronics  Research  Centre,  Department  of  Electronics  and  Electrical 
Engineering,  University  of  Glasgow,  Glasgow,G}2-8QQ,  UK,  In  the  small  size  current  and  future 
generation  semiconductor  devices  the  discrete  nature  of  the  electric  charge  emerge.  There  are  only  a  few 
hundred  impurity  atoms  and  approximately  the  same  number  of  carriers  in  a  typical  0.1  x  0.1  x  0.1  pm’ 
device.  The  random  distribution  and  statistical  fluctuation  in  the  number  of  impurity  atoms  can  produce 
variations  in  device  characteristics  that  hamper  large-scale  integration.  Variations  in  the  subthreshold 
characteristics  and  threshold  voltages  in  sub-micron  MOSFETs  caused  by  random  fluctuations  in  the  local 
impurity  distributions  have  been  clearly  demonstrated  in  simple  3D  drift-diffusion  simulations  [1].  The 
trapping  and  detrapping  of  individual  carriers  on  randomly  distributed  defect  states,  whose  pattern  is 
unique  in  any  device,  gives  rise  to  random  telegraph  noise.  Individual  carriers  can  induce  modulations  in 
the  potential  distribution  that  can  significantly  affect  the  blocking  characteristics  of  the  barrier  at  the 
source-channel  junction  of  FETs.  Furthermore,  averaged  scattering  rates  are  simply  inappropriate  in 
devices  with  countable  number  of  impurities  and  carriers,  and  the  impurity  and  carrier-carrier  scattering 
must  be  treated  as  interaction  between  individual  charges  [2].  All  these  effects  indicate  that  the  detailed 
distribution  of  the  impurities  and  the  individual  motion  of  the  carriers  must  be  considered  in  the  simulation 
of  the  semiconductor  devices  below  the  0.1  pm  limit.  A  correct  treatment  of  discrete  charges  and  their 
interaction  in  a  simulation  requires  3D  solution  of  Poisson's  equation  with  fine-grain  discretization.  This  in 
turn  requires  large  computational  resources.  A  suitable  method  to  treat  the  carrier  dynamics  must  therefore 
be  developed  to  reduce  the  computational  burden  as  much  as  possible.  The  need  to  consider  the  discrete 
nature  of  carriers  favours  methods  which  follow  the  motion  of  individual  particles  over  those  which  rely  on 
the  solution  of  partial  differential  equations.  Here  we  report  on  a  new  simple  approach  to  the  particle 
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device  simulation  based  on  Brownian  dynamics  and  the  Langevin  equation.  It  satisfies  the  above  needs 
where  the  fullblown  Monte  Carlo  method  is  not  necessary  or  too  expensive.  Apart  from  its  simplicity  our 
method  may  be  particularly  useful  for  isolating  and  studding  some  specific  aspects  of  the  carrier  dynamics 
related  to  the  "atomistic*'  nature  of  carriers  and  impurities.  It  is  more  transparent  than  the  Monte  Carlo 
technique  in  distinguishing  the  effects  due  to  the  discrete  nature  of  particles  and  their  interaction  with 
impurities,  because  these  effects  are  not  obscured  by  structure  arising  from  the  scattering  rates.  It  can  be 
also  used  for  ab-initio  noise  simulations.  First  we  will  discuss  the  relation  between  our  new  method  and 
those  currently  in  use  for  simulation  of  semiconductor  device  including  drift-diffusion,  hydrodynamic, 
Monte  Carlo  and  Cellular  Automaton.  This  will  be  followed  by  a  formal  description  of  the  Langevin 
theory  and  its  discrete  time  approximation.  The  truncation  error  and  the  convergence  of  the  discrete  time 
approximation  are  studied  in  details.  The  practical  implementation  of  the  Brownian  dynamics  approach 
will  be  described  further.  Two  modifications  of  the  Brownian  simulation  technique  have  been  developed, 
based  on  randomly  fluctuating  velocity  and  random  displacement.  Several  possibilities  for  the  choice  of 
the  magnitude  of  the  Wiener  process,  describing  both  the  velocity  and  the  displacement  fluctuation,  which 
lead  to  a  correct  Einstein  relation  will  be  presented.  Special  attention  is  also  paid  to  the  realistic  modeling 
of  the  ohmic  contacts.  First  we  have  tested  the  Brownian  approach  in  a  bulk  material  simulation  where  it 
predicts  correctly  the  carrier  temperature  and  the  diffusion  coefficient.  The  approach  has  also  been  tested 
in  a  simulation  of  a  p-n  junction  diode.  The  results  of  this  test  are  in  excellent  agreement  both  with  the 
ideal  diode  equation  and  with  the  drift-diffusion  results  obtained  from  the  commercial  simulator  MEDICI 
(Fig.  1).  The  potential  across  the  diode  obtained  from  the  Brownian  simulation  at  different  applied 
voltages  is  given  in  Fig.  2.  The  approach  will  be  further  applied  for  simulation  of  a  simplified  MESFET 
with  random  distribution  of  impurities  (Fig.3).  The  potential  distribution  in  the  middle  of  such  device  when 
both  the  carriers  and  the  impurities  are  "atomistically"  treated  is  given  in  Fig.  4. 
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TA4.  Non-Local  Impact  Ionization  Effects  of  Hot  Carriers  in  Polar  Semiconductor  Devices,  P.  Lugli, 
C.  Cianci,  A.  Di  Carlo,  and  D.  Meglio,  Dipartimento  di  Ingegneria  Elettronica,  Universitd  di  Roma  **Tor 
Vergata'\  00133  Roma,  Italy,  We  have  performed  a  Monte  Carlo  study  of  GaAs/GaAlAs  and  InP/InGaAs 
Heterojunction  Bipolar  Transistors  under  near-breakdown  conditions.  The  good  agreement  of  the 
calculated  values  of  the  multiplication  factor  with  available  experimental  measurements  allows  us  to  point 
out  the  fundamental  role  of  the  dead  space  effects  in  the  impact  ionization  phenomena  occurring  in  the 
collector.  A  first  consequence  of  our  finding  is  the  impossibility  to  extract  physical  ionization  coefficients 
from  the  measured  multiplication  factors  if  such  non-local  effects  connected  to  the  dead  space  are  not 
accounted  for.  Thus,  all  experimentally-deduced  ionization  coefficients  present  in  the  literature  for  polar 
materials  do  not  have  a  general  character  as  they  have  been  derived  from  essentially  local  models.  Such 
coefficients  are  intrinsically  dependent  on  the  specific  structure  used  for  the  measurement.  As  a  matter  of 
example,  we  will  show  that  the  use  of  the  generally  accepted  ionization  coefficients  measured  by  Bulmann, 
et  al.,  in  GaAs  photodiodes  does  not  lead  to  the  measured  multiplication  factor  of  GaAs  HBTs.  We  have, 
therefore,  developed  a  non-local  “delay”  model  which  includes  all  basic  physical  features  of  the  ionization 
process  as  extracted  from  Monte  Carlo  simulations.  The  delay  model  has  been  implemented  in  a  standard 
drift-diffusion  approach,  thus  providing  a  fast  and  reliable  device  modelling  tool  which  produces  very  well 
the  HBT  experimental  data.  Our  model  has  allowed  us  to  design  a  general  procedure  for  the  extraction  of 
physical  ionization  coefficients  from  the  multiplication  factor  of  a  variety  of  different  structures. 


P35.  Mionte  Carlo  and  Hydrodynamical  Simulation  of  a  One  Dimensional  n* — n — n*  Silicon  Diod^  O, 

Muscato\  M.V,  Fischetti,  and  R,M,  Pidatella\  ^Dipartimento  di  Matematica,Viale  Andrea  Doria,  95125 
Catania,  Italy,  ^I,B,M,  Research  Division,  T,J,  Watson  Reseach  Center,  P,0.  Box  218,  Yorktown  Heights, 
New  York  10598.  Hydrodynamical  models  are  currently  used  in  simulating  charge  carrier  transport  in 
semiconductor  devices  in  order  to  describe  high-field  phenomena  such  as  hot  electrons,  impact  ionization, 
etc.  Hydrodynamical  models  are  derived  from  the  moment  equations  of  the  Boltzmann  transport  equation 
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(BTE)  by  making  suitable  assumptions  on  the  closure  of  the  infinite  hierarchy  and  by  modeling  the  RHS's 
of  the  moment  equations  (i.e.  the  production  terms).  The  most  popular  hydrodynamical  model  is  the  one 
introduced  by  Bloetekjaer  [1]  and  subsequently  widely  investigated  by  Baccarani  and  Wordeman  [2],  in 
which  the  closure  is  achieved  by  using  the  Fourier  law  with  the  Wiedemann-Franz  expression  for  the  heat 
conductivity:  such  a  model  leads  to  serious  difficulties  when  compared  with  results  arising  from  a  direct 
numerical  simulation  of  the  BTE  [3].  Recently  Anile,  Muscato  and  Pennisi  [4], [5]  proposed  an  improved 
hydrodynamical  model  in  which  the  closure  is  achieved  by  exploiting  the  entropy  principle  [6],  and  the 
production  terms  are  modeled  as  relaxation  terms  consistently  with  the  Onsager  Reciprocity  Principle  [7]. 
In  order  to  test  such  a  model,  we  simulate  the  stationary  electron  flow  in  the  n—n — n  submicron  diode, 
which  mimics  the  channel  in  a  MOSFET  device  [8].  The  system  is  discretized  by  using  finite  differences 
and  the  1-D  box  method.  The  solution  of  the  resulting  non-linear  system  is  obtained  by  Newton’s  method 
with  the  Bank  and  Rose  [9]  damping.  In  such  a  model  the  relaxation  times  for  energy,  momentum,  energy 
flow,  shear  (appearing  in  the  production  terms)  are  determined  by  using  the  Monte  Carlo  code  Damocles^ 
[JO]:  we  also  obtained  the  macroscopic  quantities  (average  energy,  velocity,  density,  etc.),  which  are 
compared  with  those  obtained  by  the  hydrodynamical  model. 


[1]  K.  Blotekjaer,  IEEE  Trans,  on  Electron  Devices,  ED-17,  38,  (1970). 

[2]  G.  Baccarani,  M.R.  Wordeman,  Solid-State  Electronics  29,  970,  (1982). 

[3]  M.A.  Stettler,  M.A.  Alam,  M.S.  Lundstrom,  IEEE  Trans,  on  Electron  Devices  ED-40,  733,  (1993). 

[4]  A.M.  Anile,  S.  Pennisi,  Physical  Review  B  46,  13186,  (1992). 

[5]  A.M.  Anile,  O.  Muscato,  to  appear  in  Physical  Review  B  (1995). 

[6]  L  Muller  and  T.  Ruggeri,  Extended  Thermodynamics,  (Springer- Verlag,  Berlin,  1993). 

[7]  L.  Onsager,  Phys.  Rev.  37, 405,  (1931). 

[8]  C.L.  Gardner,  J.W.  Jerome,  D.J.  Rose,  IEEE  Comp.-Aided  Design  8,  501,  (1989). 

[9]  R.E.  Bank,  D.J.  Rose,  Num.  Math  37,  279  (1981). 

[10]  M.V.  Fischetti,  S.  Laux,  Phys.  Rev.  B  48,  2244,  (1993). 


P36.  A  New  Self-Consistent  2D  Device  Simulator  Based  on  Deterministic  Solution  of  the  Boltzmann, 
Poisson  and  Hole-Continuity  Equations,  Wenchao  Liang,  Neil  Goldsman  and  Isaak  Mayergoyz, 
University  of  Maryland,  Dept,  of  Electrical  Engineering,  College  Park,  MD  20742,  We  report  on  a  new 
2D  simulation  tool  which,  to  our  knowledge,  for  the  first  time  performs  actual  2D  device  modeling  by 
deterministic,  self-consistent  solution  of  the  Boltzmann  Transport  Equation  for  electrons,  the  Hole-Current 
Continuity  Equation  and  the  Poisson  Equation.  The  method  has  the  following  additional  attributes:  (i)  It 
employs  the  first  2  bands  of  the  transport  model  of  [1],  and  thereby  is  virtually  equivalent,  from  a  physical 
point  of  view,  to  performing  2-D  self-consistent  spherical  band  MC  simulations,  (ii)  The  solution  directly 
gives  the  electron  distribution  function,  electrostatic  potential,  and  the  hole  concentration  for  the  entire  2-D 
device.  Average  quantities  such  as  electron  concentration  and  electron  temperature  are  obtained  directly 
from  the  integration  of  the  distribution  function.  Impact  ionization  and  substrate  current,  as  well  as  I-V 
characteristics,  are  also  calculated  directly  from  the  distribution  function.  (iii)The  method  has  the 
advantage  of  being  significantly  faster  than  MC  simulations  and  does  not  have  statistical  noise.  The  method 
provides  more  information  than  the  hydrodynamic  approach,  and  does  not  use  empirical  mobility  models, 
(iv)  The  self-consistent,  2D  BTE  simulator  has  been  adapted  to  model  both  2D  conventional  MOSFETs 
and  SOI  MOSFETs.  Tractable  Formulation  of  the  BTE:  We  use  the  generalized  spherical  harmonic 
approach  to  transform  the  BTE  from  a  5-dimensional  differential-integral  equation  to  a  3-dimensional 
system  of  differential-difference  equations  which  is  tractable  for  numerical  evaluation[2].  With  this 
approach,  the  momentum  distribution  function  is  expressed  in  terms  of  an  infinite  series  of  spherical 
harmonics.  By  taking  advantage  of  the  recurrence  and  orthogonal  relationships  between  spherical 
harmonics,  we  derive  a  system  of  equations  for  the  expansion  coefficients  to  arbitrarily  high  order.  This 
system  is  then  truncated  after  the  4th  spherical  harmonic,  and  solved  numerically.  Numerical  Approach: 
We  developed  a  Scharfetter-Gummel  (SG)  type  discretization  for  the  BTE,  and  use  the  standard  SG 
discretization  on  the  hole-continuity  equation.  Each  equation  in  the  discretized  system  is  thereby  well- 
conditioned.  We  overcome  memory  constraints  of  the  3D  BTE  (2D  real  space,  ID  energy  space)  by  solving 
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the  discretized  BTE  with  a  fixed  point  iteration  method.  The  overall  system  is  solved  using  a  Gummel  type 
decoupled  method. 
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P37.  Self-Consistent  Solution  of  the  Poisson  and  the  Multi  Band  Boltzmann  Equations  in  Silicon, 

Surinder  P,  Singh^  Neil  Goldsman,  and  Isaak  D.  Mayergoyz,  University  of  Maryland,  Dept  of  Electrical 
Engineering,  College  Park,  MD  20742.  The  spherical  harmonic  approach  for  solving  the  BTE,  which 
represents  a  compromise  method  in  between  hydro-dynamic  modeling  and  full  zone  Monte  Carlo  (MC) 
simulation,  is  gaining  more  acceptance  [1,2].  In  this  work  we  have  developed  a  new  numerical  technique 
to  self-consistently  solve  the  multi-band  BTE  and  Poisson  equation  and  thereby  obtain  the  distribution 
function  for  energies  greater  than  3eV.  The  multi-band  band  structure  is  from  [3],  which  was  developed 
for  computationally  efficient  multi-band  MC.  The  novel  aspects  of  the  numerical  solution  include  a 
Scharfetter-Gummel  type  discretization  on  a  boundary  fitted  coordinates  (BFC)  grid  such  as  those  used  in 
computational  fluid  dynamics  (CFD)  [4].  The  previous  solutions  [2]  have  used  an  orthogonal  grid,  but  the 
numerical  solution  is  more  amenable  on  a  BFC  grid,  which,  to  the  authors  knowledge,  is  being  used  for  the 
first  time  in  semiconductor  device  simulation.  The  analytical  model  consists  of  the  BTE  for  each  band  k  (k 
=  1,2, 3, 4),  and  Poisson  equation.  The  dimensionality  of  the  BTE  is  reduced  by  projecting  it  on  a  spherical 
harmonics  basis  in  k  space,  f(r,k)  =  ^  where,  k  =  I  1,0  is  the  polar  angle,  and  0  is 


the  azimuthal  angle.  By  substitution  the  expansion  into  the  BTE,  and  applying  the  orthogonality  of 
and  changing  variables  to  F^(x,H)=  J^(x,H  +  q(l))  we  get  [2]: 
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=  0  (1).  Four  coupled  BTE  of  this  form  can  be  written  using  the 


band  structure  from  [3].  Reduction  in  computational  effort  is  achieved  by  combining  two  BTE’s  for  the 
lower  bands  into  one,  indexed  by  (/),  and  the  other  two  BTE’s  for  the  upper  bands  into  another,  indexed  by 
[//].  We  have  to  solve  the  Poisson  equation  along  with  the  BTE,  eq.  (1).  The  BTE  is  defined  on  a 
curvilinear  region  bounded  by  real  space  x  =  0,  and  H  Since  the 


boundaries  of  the  solution  domain  are  not  orthogonal  we  resort  to  a  BFC  grid  which  conforms  to  this 
boundary  and  perform  interpolations  when  writing  the  finite  difference  approximations  of  the  derivatives 
[4].  There  are  three  main  advantages  of  BFC,  (i)  the  grid  can  be  made  selectively  dense  in  B ,  (ii) 
boundary  condition  at  band  maxima  is  directly  specifiable  (iii)  during  self  consistent  calculations  the  grid 
remains  the  same.  The  equations  are  discretized  by  finite  differences,  and  using  Scharfetter-Gummel  type 
discretization  of  equation  (1).  The  individual  equations  are  solved  using  a  fixed  point  SOR  approach  while 
the  overall  system  is  solved  self  consistently  with  a  decoupled  Gummel-type  scheme.  The  approach  was 
tested  on  a  /i"  -  n  -  n'"  device  structure. 
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[3]  R.  Brunetti  and  C.  Jacoboni,  F.  Venturi,  E.  Sangiorgi,  and  B.  Ricco,  Solid-State  Electron.,  32(12),  1663, 
1989. 
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P38.  Hot-Carrier  Reliability  of  MOSFET  Isolation  Oxides  Using  3-D  Hydrodynamic  Simulation, 

Daniel  C.  Kerr,  Wenchao  Liang,  Isaak  D.  Mayergoyz,  and  Neil  Goldsman,  Department  of  Electrical 
Engineering,  University  of  Maryland,  College  Park,  MD  20742.  Hot-carrier  reliability  is  a  major  factor  in 
the  design  of  advanced  CMOS  processes.  The  isolation  oxide  is  subject  to  mechanical  stress  and  field 
boron  penetration,  which  makes  it  more  susceptible  to  hot-carrier  damage  than  the  high-quality  gate  oxide. 
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Greater  oxide  damage  at  the  channel  edges  causes  the  current  path  to  redistribute,  influencing  the  hot- 
carrier  lifetime.  Experimental  studies  have  difficulty  resolving  the  location  of  oxide  degradation.  Estimates 
of  hot-carrier  lifetime  from  substrate  current  or  threshold  voltage  shift  measurements  may  be  inaccurate 
due  to  the  3-D  redistribution  of  the  current  path,  which  changes  the  rate  of  degradation.  Existing  numerical 
studies  have  used  2-D  simulation.  This  paper  announces  a  new,  3-D  hydrodynamic  T-CAD  tool  which  can 
predict  the  detailed  accumulation  of  interface  and  oxide  traps  over  time  and  their  influence  on  the  hot- 
carrier  lifetime.  SIMaster  solves  the  hydrodynamic  model  equations  formulated  in  terms  of  new  state 
variables  especially  tailored  for  the  equation-  and  space-decoupled  fixed-point  iteration  technique[l,2]. 
The  gate  current,  interface  state  generation,  and  oxide  trap  generation  are  computed  based  on  an  electron 
distribution  function  fit  to  the  average  electron  energy  at  successive  time  steps  (Fig.  1).  This  model  has 
been  successfully  used  to  model  MOSFET  reliability  as  well  as  EPROM  and  FLASH  cell  programming 
[3].  The  model  has  been  applied  to  a  0.25  ^rn  LOCOS-isolated  MOSFET  (Fig.  2).  The  stress  condition 
wasV^  =  4.5V,  Vq  =  5.5V  for  a  stress  time  T^  =  10*s.  The  electron  temperature  in  a  2-D  slice  at  the  center  of 
the  channel  is  shown  in  Fig.  3.  The  electron  temperature  in  a  slice  immediately  beneath  the  gate  oxide  is 
shown  in  Fig.  4.  Fig.  5  shows  the  doping  concentration,  electron  concentration,  and  electron  temperature 
in  the  channel.  The  computed  interface  and  oxide  traps  for  each  time  step  are  shown  in  Fig.  6,  which  is  a 
line  across  the  center  of  the  channel. 
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[3]  J.  Z.  Peng,  Q.  Lin,  P.  Fang,  M.  Kwan,  S.  Longcor,  and  J.  Lien.  "Accurate  Simulation  of  EPROM  Hot- 
Carrier  Induced  Degradation  using  Physics-Based  Interface  and  Oxide  Charge  Generation  Models," 
International  Reliability  Physics  Symposium,  pp.  154-160,  1994. 


P39.  3-D  Device  Simulation  Using  Intelligent  Solution  Method  Control,  Daniel  C.  Kerr  and  Isaak  D. 
Mayergoyz,  Department  of  Electrical  Engineering,  University  of  Maryland,  College  Park,  MD  20742. 
Numerical  simulation  of  semiconductor  devices  requires  the  solution  of  a  large,  nonlinear,  coupled  system 
of  discrete  equations.  The  discretized  variables  of  the  model  equations  (Poisson,  electron  and  hole  current- 
continuity,  etc.)  are  coupled  both  in  equation  space  and  in  real  space  (see  Fig.  1).  Several  solution  methods 
which  treat  the  coupling  differently  have  been  extensively  developed:  full-coupled  (Newton),  decoupled  in 
equation  space  (Gummel),  and  decoupled  both  in  equation  and  real  space  (fixed-point  or  Jacobi  iteration) 
(see  Fig.  1(a)).  Each  solution  method  can  be  explained  by  using  Fig.  1(b).  The  Newton  method  includes  all 
the  coupling  in  equation  and  real  space.  The  Newton  method  gathers  together  the  variables  from  all  the 
equations  at  all  mesh  points  into  one  large  matrix  equation  and  solves  it  iteratively.  The  Gummel  method 
gathers  up  the  variables  at  all  mesh  points  for  each  equation  into  different  matrix  equations  and  solves  them 
iteratively.  The  fixed-point  method  doesn't  use  matrices  at  all.  Instead,  a  single  equation  is  solved  and 
iterations  are  performed  over  all  mesh  points  and  all  equations.  Finally,  the  point-Newton  method  takes 
the  variables  for  all  the  equations  at  only  one  mesh  point  into  a  small  matrix  and  iterates  over  all  mesh 
points.  Each  of  these  methods  has  its  strengths  and  weaknesses  for  different  simulation  problems.  For 
Newton's  method,  the  advantage  is  a  quadratic  rate  of  convergence;  disadvantages  are  it  requires  a  good 
initial  guess  (not  globally  convergent)  and  a  large  computer  memory.  Gummel's  method  reduces  the 
computer  memory  requirements,  but  this  method  will  slow  or  fail  for  tightly  coupled  equations  (high  bias 
conditions).  The  fixed-point  iteration  method  is  globally  convergent,  not  memory  intensive,  and  inherently 
parallel;  however,  it  suffers  the  problems  of  Gummel's  method  and  slow  convergence  when  the  mesh 
points  are  tightly  coupled.  The  new  solution  algorithm  is  illustrated  in  Fig.  2.  After  an  initial  guess  is 
calculated,  the  degree  of  equation  and  space  coupling  between  discretized  variables  is  calculated,  and  the 
mesh  points  are  partitioned  into  groups.  The  best  available  solver  is  applied  to  each  group,  and  the 
iteration  is  over  the  groups.  This  approach  avoids  the  weaknesses  of  each  solver  by  matching  its  strengths 
to  a  narrowly  tailored  problem.  The  degree  of  equation  coupling  is  measured  by  deviation  from  charge 
neutrality  or  inversion.  For  example,  in  a  MOSFET,  this  assigns  the  channel  region  to  equation-space 
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coupled  solvers.  The  substrate  points  are  quasi-neutral,  so  they  are  assigned  to  equation-space  decoupled 
solvers.  The  degree  of  space  coupling  is  measured  by  mesh  nonuniformity  and  by  strong  current  flows. 
Thus,  the  channel  region  will  be  solved  using  Newton’s  method,  the  depletion  regions  surrounding  the 
source  and  drain  will  be  solved  using  Gummel's  method,  and  the  rest  by  the  fixed-point  iteration  method. 
This  efficient  and  intelligent  partitioning  saves  CPU  time  compared  to  any  single  method  alone.  This 
selective  coupling  approach  has  been  integrated  into  the  general  3-D  device  simulation  program 
SIMAsTER  [1].  Newton’s  method  can  be  applied  to  very  large  meshes  on  engineering  workstations  by 
restricting  it  to  the  critical  coupled  points,  which  speed  improvements  over  other  methods.  Rates  of 
convergence  and  solution  times  will  be  given,  and  compared  to  standard  simulation  programs  (PISCES-II). 


Supported  by  the  Semiconductor  Research  Corporation  contract  SJ-377. 
[1]  D.C.  Kerr  and  1.  D.  Mayergoyz,  to  appear,  1995. 


P40.  3D  Parallel  Monte  Carlo  Simulation  of  GaAs  MESFETs,  S.  Pennathur,  Can  K.  Sandalci,  Cetin  K. 
KoCy  and  S.M.  Goodnick,  Department  of  Electrical  and  Computer  Engineering,  Oregon  State  University, 
Corvallis,  OR  9733 L  We  have  investigated  three-dimensional  (3D)  effects  in  sub-micron  GaAs  MESFETs 
using  a  parallel  Monte  Carlo  device  simulator,  PMC-3D  [1].  The  parallel  algorithm  couples  a  standard 
Monte  Carlo  particle  simulator  for  the  Boltzmann  equation  with  a  3D  Poisson  solver  using  spatial 
decomposition  of  the  device  domain  onto  separate  processors.  The  3D  Poisson  equation  is  solved 
iteratively  in  parallel  using  a  red-black  coloring  scheme  for  up  to  one  million  grid  points  on  distributed 
memory  multiprocessor  machines.  For  realistic  3D  device  structures,  we  find  that  the  main  performance 
bottleneck  is  the  Poisson  solver  rather  than  the  Monte  Carlo  particle  simulator  for  the  parallel  successive 
overrelaxation  (SOR)  scheme  employed  in  [1].  Thus  parallel  multigrid  algorithms  are  studied  and 
compared  to  the  previous  SOR  implementation.  The  frequency  dependent  small  signal  parameters  of  a 
MESFET  device  have  been  studied  as  a  test  case  of  the  3D  algorithm.  Parameters  such  as  the 
transconductance  and  output  impedance  are  calculated  through  frequency  domain  analysis  of  the  transient 
current  response  to  small  electrode  voltage  steps.  The  scaling  properties  of  the  small  signal  parameters 
have  been  simulated  for  both  the  gate  width  in  the  third  dimension  as  well  as  the  gate  length.  Deviations 
from  the  ideal  linear  scaling  of  the  device  gain  with  lateral  scaling  of  the  gate  width  are  found  associated 
with  the  3D  field  and  particle  distribution  in  the  structure. 


[1]  U.A.  Ranawake,  C.  Huster,  P.M.  Lenders,  and  S.M.  Goodnick,  IEEE  Trans,  on  CAD,  Vol.  13,  No.  6, 
712(1994). 


P  41.  Hydrodynamic  Simulations:  The  Forced  Notched  Oscillator,  and  Channel  Geometry  in  the 
MESFET,  Gui-Qiang  Chen  and  Joseph  Jerome\  Department  of  Mathematics,  Northwestern  University, 
Evanston,  IL  60208;  Chi-Wang  Shu,  Division  of  Applied  Mathematics,  Brown  University,  Providence,  RI 
02912.  We  introduce  a  novel  two  carrier  hydrodynamic  model,  which  has  the  capability  of  incorporating 
higher  dimensional  geometric  effects  into  a  one  dimensional  model,  if  symmetry  is  present.  This  abstract 
presents  two  special  cases  of  the  model,  the  first  of  which  will  be  illustrated  on  the  following  two  pages 
through  simulations.  For  the  GaAs  device  in  the  notched  oscillator  circuit,  considered  by  the  authors  of 
[3],  the  model  reduces  to  that  first  considered  by  Blotekjaer  in  [1],  as  a  two  valley  hydrodynamic  moment 
model.  However,  the  critical  form  of  the  couplings  is  not  discussed  in  [1],  and  is  studied  in  detail  here. 
We  do  not  employ  the  actual  circuit,  but  use  an  analogue  sinusoidal  forcing  term  for  the  bias  at  the  drain, 
closely  correlated  to  the  resonant  frequency  of  the  circuit,  which  can  be  computed  by  analytical  methods. 
We  find  that  our  concentration  curves,  illustrated  for  both  valleys  (l=r,  2=L),  share  the  structural  features 
of  [3],  where  the  curves  were  computed  by  Monte-Carlo  simulation.  It  is  also  apparent  that  the  coupling 
threshold  in  the  energy  equation  for  the  F  valley  is  crucial  for  the  inter-valley  transfer  between  the 
carriers.  We  have  performed  simulations  for  different  coupling  thresholds  with  different  results  (not  shown 
in  the  abstract).  In  fact,  it  is  simply  inadequate  information  about  the  couplings  which  limits  the  valleys  to 
two  in  this  modeling.  We  acknowledge  with  gratitude  the  extremely  helpful  advice  of  Umberto  Ravaioli 


40 


for  this  part  of  the  work.  If  the  one  carrier  MESFET  channel  is  geometrically  idealized  by  a  nozzle 
configuration,  the  model  reduces  to  a  one  dimensional  perturbed  hydrodynamic  model  (one  carrier).  The 
perturbation  is  a  geometric  source  term,  proportional  to  the  derivative  of  the  cross  sectional  area  of  the 
channel.  It  has  been  found  in  transonic  nozzle  flow  that  this  term  assumes  an  important  role  in  predicting 
possible  shock  structure,  and  we  study  its  role  in  semiconductor  modeling.  Our  algorithms  are  extensions 
of  those  described  in  [2]. 


[1]  K.  Blotekjaer,  Transport  equations  for  electrons  in  two-valley  semiconductors,  IEEE  Trans.  Electron 
Devices  ED47  (1970)  38-47. 

[2]  J.W.  Jerome  and  C.-W.  Shu,  Energy  models  for  one-carrier  transport  in  semiconductor  devices,  in  IMA 
Volumes  in  Mathematics  and  Its  Applications.  Vol.59,  W.  Coughran,  J.  Cole,  P.  Lloyd  and  I.  White,  eds., 
Springer-Verlag  1994  pp.  185-207. 
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P42.  Transient  Analysis  of  Silicon  Devices  Using  the  Hydrodynamic  Model,  Luigi  Colalongo,  Marina 
Valdinoci,  Antonio  Gnudi,  and  Massimo  Rudan,  Dipartimento  di  Elettronica,  Universitd,  di  Bologna,  Viale 
Risorgimento  2,  40136  Bologna,  Italy,  The  analysis  of  the  switching  behaviour  of  submicron  devices 
brings  about  the  necessity  of  extending  the  solution  of  the  Hydrodynamic  Model  to  the  transient  case.  The 
implementation  of  such  model  has  been  carried  out  and  a  few  examples  of  simulation  are  presented  here, 
showing  the  velocity-overshoot  of  a  ballistic  diode  and  the  temperature  spread  in  the  drain  region  of  a 
realistic  MOS  device.  The  Hydrodynamic  Model  (HD)  in  the  transient  case  reads  [1]: 

—<ii\(e,gradcp)  =  q(p-n  +  N*-  N^), 

dnldt-dhliJJq)=-U, 

dp! dt  +div{J^I q)=U,  (1) 

d{n(0„)ldt+div\-k^gradT^-{5l2)kgTJJq'\=B*J^  -(oJJ-  n{a„  - 

(2) 

d(,p(Op )  /dt  +  div\-k^gradTp  +  (5 /  2)kJ^  J/q]=E»J^- 0)^U  -  p(o)^  -  , 

(3) 

where  J.  =  qD,gradn+^vngrad(k,T„/q-(J5),  =  -qD^gmdp  +  qp^gradik^T^/ q  -(p),  T,  are 

the  temperature,  momentum-relaxation  time,  and  energy-relaxation  time  of  the  carriers,  and 
(0=(XI2)mv^  +Ol2)kJ  is  the  average  kinetic  energy.  The  remaining  symbols  have  the  usual 
meaning.  This  set  of  equations  is  discretized  in  space  using  the  so-called  Box  Integration  Method,  The  time 
derivatives  are  discretized  using  the  Gear  method  which  is  an  A-stable,  L-stable  method  for  stiff  differential 
equations.  The  discretization  yields  at  each  time  step  a  system  of  5N  algebraic  equations,  where  N  is  the 
size  of  the  discretization  grid  [1].  The  system  of  non-linear  equations  is  solved  using  the  following  scheme: 
first,  the  system  of  Poisson,  electron-  and  hole-continuity  equations  is  solved  by  a  coupled  Newton  method, 
then  the  temperature  is  updated  solving  the  energy-balance  equations  for  electrons  and  holes.  The  external 
loop  is  repeated  until  the  global  convergence  is  reached,  and  the  whole  procedure  is  repeated  at  each  time 
step.  The  solution  scheme  depicted  above  is  easier  to  implement  than  the  full  Newton  method  coupling  all 
equations;  it  has  proven  very  stable  and,  supplemented  with  a  quadratic  projection  algorithm,  rather  fast  as 
well.  The  term  E  •  J  is  treated  by  the  vectorial  identity  E  =  ’-div{(pj)  +  (pdiv{J)[3L  this  leading 

to  This  expression  is  easy  to  handle  as  the  time  dependence  is  not  explicit: 

therefore,  in  the  transient  case  there  is  no  need  to  account  for  the  time  derivatives  in  div(J)  and  div(^  J). 
The  simulations  have  been  carried  out  using  a  two-dimensional  version  of  the  device-analysis  program 
HFIELDS  [1],  supplemented  with  the  method  described  above,  first  on  a  0.5  pm,  n-n-n  ballistic  diode 
having  a  0.35  pm  n-region.  The  diode  is  initially  at  equilibrium,  then  it  is  biased  by  a  linear  ramp  which 
brings  the  voltage  drop  between  the  contacts  from  0  to  2  V  in  1  psec;  the  voltage  on  this  contact  remains  at 
2  V  for  1  psec  after  the  end  of  the  ramp.  It  is  worth  noting  that  the  duration  of  the  ramp  is  about  one  order 
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of  magnitude  larger  than  the  carrier  momentum-relaxation  time;  for  this  reason  the  terms  dj/di  are 
dropped  from  the  transport  equations.  On  the  contrary,  the  time  derivatives  are  retained  in  (2,3)  because  the 
carrier  energy-relaxation  times  are  larger  [2].  In  any  case  the  ramp  is  slow  enough  to  allow  for  the  steady- 
state  approximation  leading  to  Poisson's  equation.  Fig.  1  shows  the  normalized  temperature  at  two 
different  times  1  psec  and  I  psec.  The  first  one  is  of  the  same  order  as  X(0  and  corresponds  to  the  instant  at 
which  the  ramp  reaches  2  V;  at  r  =  1  psec  the  transient  is  extinguished.  The  figure  exhibits  the  delay  in  the 
temperature  growth  due  to  the  time  derivatives  in  (2,3).  In  Fig.  2  a  velocity  overshoot  is  observed  at  r  =  1 
psec,  due  to  the  fact  that  the  momentum-relaxation  time  is  smaller  by  about  one  order  of  magnitude  than 
the  energy-relaxation  time.  Because  of  this,  the  velocity  increase  is  considerably  faster  than  that  of  energy, 
as  is  seen  in  Figs.  1  and  2.  As  the  carrier  temperature  grows,  the  velocity  overshoot  becomes  less 
important.  The  features  of  the  momentum  and  energy-relaxation  times  are  reflected  into  the  mobility  model 
used  in  the  HD  equations,  /i(fi))  =  +0^(6)  “■  intervals 

(t  <  T^),  while  the  temperature  is  still  growing  the  mobility  is  in  fact  larger  than  in  steady-state 
conditions.  Figs.  3  and  4  show  the  electron  concentration  and  the  electric  field,  respectively,  at  the  same 
time  intervals  and  at  equilibrium.  Next,  more  simulations  have  been  carried  out  on  a  realistic  n-type,  0.8 
pm  channel  MOS  transistor.  The  source  and  gate  voltages  were  set  at  0  and  1  V,  respectively.  Starting  from 
equilibrium,  the  drain  voltage  was  first  brought  from  0  to  3  V  using  a  1-psec  linear  ramp,  then  kept  at  3  V 
for  1  psec.  Figs.  5  and  6  show  the  normalized  temperature  at  ?  =  1  psec,  i.e.,  when  the  ramp  reaches  3  V 
and,  respectively,  at  the  end  of  the  transient.  In  the  direction  of  the  channel  the  results  are  qualitatively 
similar  to  those  of  the  previous  case;  the  temperature  peak  is  initially  localized  in  the  channel  region  while, 
after  the  end  of  the  ramp,  it  starts  spreading  to  eventually  reach  the  steady-state  profile  of  Fig.  6.  In 
addition,  thanks  to  the  two-dimensional  analysis,  the  spread  of  the  energy  from  the  surface  to  the  bulk  is 
also  evident  in  the  figure. 
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P43.  Numerical  Simulation  of  Thermal  Behavior  in  a  Trench  Emitter  Insulated  Gate  Bipolar 
Transistor  (IGBT),  LSabesan,  P.  Mawby,  M.S.  Towers,  K.  Board,  P.  Waind"^,  Department  of  Electrical  & 
Electronic  Engineering,  University  of  Wales  Swansea,  Singleton  Park,  Swansea  SA2  8PP,  UK,;  *  GEC~ 
Plessey  Semiconductors,  Carhome  Road,  Lincoln,  LNl  ISG,  UK,  The  modelling  of  thermal  effects  in  an 
IGBT  structure  is  presented.  The  semiconductor  equations  are  solved  in  two  dimensions  with  physical 
effects  such  as  carrier-carrier  scattering  mobility  and  SRH  and  Auger  recombination  included.  I-V 
characteristics,  carrier  concentrations  and  temperature  profiles  have  been  investigated.  Among  the  various 
power  Bipolar  and  MOS  transistors  the  IGBT  provides  better  performance  with  respect  to  low  on- 
resistance  and  low  gate  drive  requirements.  Further,  IGBTs  offer  better  high-temperature  forward 
conduction  characteristics  by  being  less  sensitive  to  temperature  variation  in  comparison  to  other  power 
devices  [2,3].  This  feature  makes  the  device  attractive  for  applications  in  which  high  ambient  temperatures 
may  be  encountered.  The  main  drawback  in  the  device  is  its  tendency  to  latchup  at  high  current  levels. 
This  limits  the  operation  of  the  device  within  a  range  determined  by  various  structural  and  material 
parameters.  In  order  to  overcome  this,  a  new  device  structure,  the  trench  emitter  configuration,  has  been 
considered  and  is  reported  here.  The  schematic  structure  of  the  device  is  shown  in  figure  1.  This  structure 
is  seen  to  be  slightly  different  from  conventional  IGBT  structures  shown  in  references  [3-8].  In  order  to 
eliminate  the  possibility  of  punch-through  in  the  device,  a  buffer  layer  has  been  introduced  between  the 
substrate  and  the  n-base  region.  The  thick  lightly  doped  epitaxial  layer  is  necessary  to  support  high 
voltages  in  the  forward  blocking  mode.  However,  it  also  contributes  to  a  large  on-state  resistance.  This 
paper  analyses  the  dependence  of  the  operating  characteristics  and  the  internal  carrier  distributions  on  the 
ambient  temperature  under  isothermal  conditions  and  considers  the  effects  of  Joule  heating  and 
recombination  in  the  non-isothermal  case.  The  simulations  were  carried  out  using  the  2-D  Swansea 
University  Device  Simulator  (SUDS).  The  resulting  spatial  current  density,  electric  field  and 
recombination  effects  in  the  device  were  used  to  calculate  the  heat  dissipation  distribution.  The  heat 
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equation  was  then  solved  to  provide  the  temperature  profile  self  consistently.  The  temperature  dependence 
of  thermal  conductivity  in  silicon  was  also  considered  in  the  calculation.  The  low  on-resistance  essential  in 
all  power  devices  is  achieved  in  IGBTs  through  high  injection  of  electrons  and  holes  into  the  n-base  drift 
region  thereby  causing  conductivity  modulation.  The  simulated  output  characteristics  are  shown  in  figure 
2,  together  with  measured  results.  As  can  be  seen  a  very  good  qualitative  and  quantitative  agreement  has 
been  obtained.  Figure  3  shows  the  influence  of  ambient  temperature  on  the  output  characteristics  at  a  gate 
bias  of  15  V  under  simple  isothermal  conditions.  This  shows  the  voltage  drop  increasing  with  temperature 
due  to  the  increase  of  resistance  from  0.01  ohm  at  300K  to  0.03  ohm  at  500K.  The  change,  however  is 
relatively  small  compared  to  that  in  MOSFET  [2].  Similar  findings  were  observed  in  measurements  and 
reported  in  references  [2], [3].  Figure  4  shows  the  effect  of  temperature  variation  on  the  distribution  of 
carrier  concentration  through  a  vertical  cut  at  x  =  10  microns  with  biases  of  15  V  on  the  gate  and  SV  on  the 
anode.  Reduction  in  electron  and  hole  concentrations  can  be  seen  with  increasing  temperature.  Any 
increase  in  temperature  severely  reduces  the  mobility  of  the  carriers,  more  significantly  in  the  drift  region. 
Figure  5  shows  the  effect  of  temperature  variation  on  the  hole  mobility  through  the  same  vertical  cut  with 
the  same  biases  as  stated  for  figure  4.  This  is  the  dominant  mechanism  for  the  reduction  of  current  at 
higher  temperatures.  The  temperature  distribution  throughout  the  domain  caused  by  self  heating  is  shown 
in  figure  6.  A  small  hot  spot  can  be  seen  near  the  gate-emitter  region.  This  is  caused  mainly  by  high  Joule 
heating  due  to  locally  increased  electric  field  and  current  density. 
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P44.  Monte  Carlo  Simulation  of  a  Submicron  MOSFET  Including  Inversion  Layer  Quantization, 
J,B,Rolddn,  F.Gdmiz,  JA.Ldpez-Villanueva  and  J,E.Carceller,  Departamento  Electronica  y  Tecnologia 
ComputadoreSy  Universidad  de  Granada^  18071  Granada,  Spain.  A  simulation  of  a  submicron  MOSFET 
by  Monte  Carlo  method  including  inversion  layer  quantization  has  been  done.  The  simulator  has  been 
applied  to  the  study  of  electron  transport  in  normal  operation  conditions.  The  device  is  divided  into  a 
convenient  number  N  of  smaller  channels  (pchannels)  of  unknown  length.  The  difference  between  the 
pseudofermi  levels  is  fixed  in  the  extremes  of  each  ^channel.  The  one-dimensional  Poisson  and 
Schrodinger  equations  are  then  self-consistenly  solved  in  each  of  these  points.  The  length  of  each  fichannel 
is  iteratively  calculated  by  imposing  the  current  continuity  along  the  MOSFET  channel.  Both  drift  and 
diffusion  components  are  taken  into  account  in  a  modified  charged  sheet  model  in  order  to  evaluate  the 
drain  current.  The  resultant  potential  distribution  is  used  as  the  starting  point  to  iteratively  solve  the  two- 
dimensional  Poisson  equation.  Following,  both  the  current  and  length  of  each  ^channel  are  again 
calculated.  This  procedure  is  repeated  until  a  convergence  criteria  is  reached.  An  adaptive  grid  is  used 
along  the  channel.  In  each  point,  the  electron  scattering  rates,  and  the  longitudinal  electric  field  are 
evaluated.  To  proceed  in  the  Monte  Carlo  simulation,  a  great  number  of  electrons  are  introduced,  one  by 
one,  in  the  channel  from  the  source.  Then,  they  move  towards  the  drain  drifted  by  the  longitudinal  electric 
field,  undergoing  phonon.  Coulomb  and  surface-roughness  scattering.  Different  electron  transport 
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properties  are  calculated  from  the  mean  time  each  electron  spends  in  every  grid  zone.  Some  results  for  a 
0.2p.m  channel  length  MOSFET  are  shown  in  the  attached  figures. 


TA5.  Analysis  of  Anomalous-Negative  Differential  Resistance  in  Diode  Breakdown  Simulation  Using 
Carrier  Temperature  Dependent  Impact  Ionization,  Edwin  C.  Kan,  Zhiping  Yu  and  Robert  W.  Dutton, 
AEL  201,  Stanford  University,  Stanford,  CA  94305.  In  recent  years,  the  drift-diffusion  (DD)  model  has 
been  extended  to  the  energy  transport  (ET)  model  and  the  hydrodynamic  (HD)  model  to  account  for  the 
elevated  average  energy  in  the  carrier  systems  at  high  fields.  Inclusion  of  electron  and  hole  temperatures  as 
state  variables  allows  transport  coefficients  such  as  mobility  and  impact  ionization  rate  as  functions  of  the 
carrier  temperatures  instead  of  the  local  electric  field  to  remove  the  thermal  equilibrium  approximation  [1]. 
However,  this  new  parametrization  of  the  transport  coefficients  has  received  controversial  criticisms. 
Although  much  success  has  been  achieved  in  fitting  specific  experiment  measurements,  especially  in 
MOSFET  substrate  current  cases  (e.g.,  [2]),  it  has  also  been  established  that  the  distribution  function,  in 
particular  the  tail  part  that  determines  the  impact  ionization  and  MOSFET  gate  current,  is  not  well 
characterized  by  the  carrier  concentration  and  average  energy  alone  [3,  4].  Nevertheless,  if  the  impact 
ionization  rates  are  chosen  as  functions  of  the  carrier  temperatures,  for  the  first  time  we  have  observed  an 
anomalous  negative  differential  resistance  (NDR)  in  a  simple  pn  junction  breakdown  simulation. 
Mechanisms  that  cause  this  NDR  are  analyzed  below.  The  reverse-bias  diode  IV  curves  using  field  and 
energy  dependencies  for  the  mobility  p  and  the  impact  ionization  rate  a  are  shown  in  Fig.  1.  The  simulation 
is  performed  by  PISCES-2ET  with  curve  tracing  techniques  [5]  to  capture  NDR.  The  DD  model  with  a(F) 
does  not  have  any  NDR  region  and  the  breakdown  voltage  is  close  to  the  measured  data.  At  the  onset  of 
breakdown,  which  is  highlighted  in  Fig.  2,  the  ET  model  with  a(TJ  and  |i(r^)  already  shows  a  small  NDR, 
while  the  ET  model  with  a(T^)  and  does  not.  This  is  due  to  the  cooling  mechanism  of  impact 
ionization  [6]  that  gives  a  positive  feedback  for  enhancement  of  pfTc),  since  asymptotically  fJT^  ^  for 


velocity  saturation.  This  small  NDR  will  disappear  if  \x(F)  is  used  or  the  cooling  mechanism  is  not 
accounted  for  in  the  carrier  energy  balance  equation.  In  addition,  significant  NDR  in  the  high  current 
region  appears  in  both  ET  models  as  long  as  both  aJTJ  and  afT^  are  used.  This  NDR  cannot  be 
attributed  to  the  previous  cooling  mechanism  since  the  feedback  is  removed  ill  the  case  of  \x(F),  but  is 
possibly  caused  by  the  Joule-heat  feedback  to  the  carrier  temperature  or  the  dislocated  afTJ  and  afT^).  As 
can  be  seen  in  Fig.  3,  in  the  depletion  region  the  peaks  of  and  T^are  offset  by  a  small  distance  [6].  The 
multiplication  factor  M  defined  by  the  ratio  of  outgoing  and  incoming  electron  (or  hole)  currents  in  the 


depletion  region  can  be  expressed  as  [7]:  ^  ^  |  o^n  ^ 


If  (X^(x)  =  (Xp(x),  the  breakdown  condition  as  M  ^  is  simply  J  =  1  and  no  NDR  is  possible. 


However,  if  OC^(x)  ^  (Xp(x),  the  reverse  breakdown  current  is  affected  by  the  convolution  of  (X^  and 
OCp,  and  NDR  can  exist  based  for  different  and  profiles.  More  details  will  be  presented  at  the 
conference.  Besides,  it  can  be  observed  that  the  breakdown  voltage  predicted  by  (X(T^)  is  much  larger 

than  that  by  a(F)[9].  This  can  be  understood  from  two  aspects.  First,  when  only  one  carrier  is  considered 
(as  in  nin  cases  of  Fig.  4),  the  generation  rate  calculated  by  a(r^)  is  smaller  than  a(F)  due  to  the  latent  effect 
from  dF/dx  (note  that  in  the  bulk  case,  (X{T^  =  OC{F)  by  definition).  Second,  if  the  peaks  of  a^  and  a^  are 

dislocated,  the  integral  in  (1)  will  be  damped  when  a^  >  a^  (similar  calculations  hold  for  holes  too). 
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TA6.  Statistical  Enhancement  of  Terminal  Current  Estimation  for  Monte  Carlo  Device  Simulation, 
P.D.  Yoder y  D.  Krumbeiriy  K,  Gartner,  and  W.  Fichtner,  Swiss  Federal  Institute  of  Technology,  Integrated 
Systems  Laboratory,  Gloriastrasse  35,  CH-8092  Zurich,  Switzerland,  To  account  for  capacitative  effects 
and  to  improve  upon  the  poor  convergence  properties  of  the  particle  counting  method,  terminal  current 
estimators  may  be  devised  which  make  use  of  information  about  all  particles  within  a  given  device.  One 
such  formalism  was  developed  by  Shockley  [1]  and  Ramo  [2],  relating  the  currents  induced  on  an  arbitrary 
number  of  terminals  to  the  motion  of  charges  in  multiple  dimensions.  For  use  in  semiconductor  device 
simulation,  this  approach  has  several  drawbacks.  We  present  theory  and  applications  of  a  new  generalized 
domain  integration  technique  for  Monte  Carlo  simulation,  which  in  contrast  1)  accelerates  convergence  by 
minimizing  the  terminal  current  estimator  variances,  2)  is  valid  for  time-dependent  or  steady-state 
calculations,  3)  differentiates  between  electron,  hole  and  displacement  currents,  and  4)  reduces  to  the 
Ramo-Shockley  theorem  as  a  special  use.  Simulated  charge  density  and  velocity  information,  such  as 
shown  in  Fig.  1,  is  integrated  against  a  special  weighting  function  throughout  the  domain  for  each  bias 
condition,  to  obtain  I-V  curves  like  the  one  shown  in  Fig.  2.  Figure  2  also  demonstrates  the  accelerated 
convergence  rates  for  source,  drain  and  substrate  currents  which  may  be  achieved  with  this  technique. 
Finally,  we  show  how  this  method  may  be  applied  to  “window”-type  Monte  Carlo  simulation. 
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TA7.  New  Approach  to  Hot  Electron  Effects  in  Si-MOSFETs  Based  on  an  Evolutionary  Algorithm 
Using  a  Monte  Carlo-Like  Mutation  Operator,  J,  Jakumeit,  U,  Ravaioli,  K.  Hess,  Beckman  Institute, 
University  of  Illinois,  Urbana,  IL  6I80L  A  new  approach  to  study  gate  and  substrate  current  in  Si- 
MOSFETs  is  introduced,  which  is  based  on  a  mixture  of  evolutionary  optimization  algorithms  and  Monte 
Carlo  technique.  The  goal  is  to  develop  an  algorithm  for  the  investigation  of  hot  electron  effects  in  Si- 
MOSFETs,  which  is  less  complex  and  computationally  expensive  than  a  full  band  Monte  Carlo  program, 
yet  gives  more  precise  information  than  lucky  electron  or  electron  temperature  models  and  can  thus  be 
used  as  a  mediator  between  precise  theory  and  experimental  results.  To  achieve  this  goal,  use  is  made  of 
an  evolutionary  optimization  algorithm  EA  [1]  to  search  for  electron  distribution  functions  x),  which 
fit  a  given  substrate  current  will  usually  be  an  experimental  result.  The  EA  searches  for  the 

correct  f[E,x)  by  starting  from  a  set  of  initial  guesses  for  f(E,x),  In  each  iteration  new  guesses  are  created 
by  modifying  the  old  guesses  with  help  of  a  Monte  Carlo-like  mutation  operator.  From  the  old  and  new 
guesses  a  new  set  of  distribution  functions  is  selected  by  simulated  annealing-like  selection  criteria.  In  this 
way  the  set  of  distribution  functions  converges  towards  an  appropriate  solution.  The  algorithm  stops,  if 
calculated  from  x)  lies  within  0.1%  of  The  results  of  this  new  approach  are  distribution 

functions  for  the  region  of  interest  in  terms  of  gate  and  substrate  current,  which  can  be  compared  to  results 
of  a  full  band  Monte  Carlo  method  with  enhancement  of  the  high  energy  tail  of  the  distribution.  If  is 

not  consistent  with  the  physical  model  inside  the  mutation  operator,  the  algorithm  does  not  converge  or 
gives  unphysical  results.  An  important  advantage  of  the  EA  is  that  even  very  simplified  mutation  operators 
will  drive  the  search  for  the  distribution  function  in  the  correct  direction.  The  use  of  a  goal  function  such  as 
x))  guarantees  a  reasonable  estimate  particulary  for  the  portion  of  the  distribution  function  that 

is  responsible  for  ,  The  energy  distribution  obtained  in  this  way  can  then  be  used  to  assess  the  validity 

of  simpler  models  (e.g.  electron  temperature)  and  can  also  be  directly  compared  to  results  of  full  band 
Monte  Carlo  simualtions.  First  results  indicate  that  this  method  can  also  be  used  to  derive  relationships 
between  gate  and  substrate  currents,  that  are  superior  to  the  ones  derived  from  electron  temperature  and 
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lucky  electron  models.  Our  new  method  may  therefore  bridge  the  gap  between  these  simplified  models  and 
full  band  approaches  for  considerations  of  hot  electron  effects  on  device  reliability. 
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WAl.  Absorbing  Boundary  Conditions  for  Quantum  Evolution  Equations,  Anton  Arnold,  Department 
of  Mathematics,  Purdue  University,  West  Lafayette,  IN.  When  numerically  solving  a  whole-space 
evolution  problem,  the  computation  has  to  be  confined  to  a  finite  domain  by  introducing  artificial  boundary 
conditions.  Numerically  transparent  or  absorbing  boundary  conditions  (ABC)  have  first  been  derived  by 
Engquist  &  Majda  for  hyperbolic  systems  by  requiring  that  outgoing  waves  can  leave  the  computational 
domain  without  being  reflected  back  in.  This  talk  will  be  concerned  with  the  construction,  analysis,  and 
discretization  of  ABC's  for  the  (transient)  Schrodinger  equation  and  the  (relaxation-time)  Wigner  equation. 
Under  the  assumption  of  a  homogeneous  medium  (i.e.  constant  potential)  outside  of  the  computational 
domain  one  can  derive  a  transparent  BC  for  the  Schrodinger  equation.  In  ID  it  relates  the  normal  space 
derivative  of  the  wave  function  at  the  boundary  to  its  fractial  (1/2)  time  derivative,  which  can  be  written  as 
a  memory  term  involving  the  past  evolution.  In  2D  the  transparent  BC  is  also  non-local  in  space.  For 
numerical  efficiency,  it  therefore  has  to  be  approximated  by  an  ABC  that  is  at  least  local  in  space.  When 
using  the  Crank-Nicolson  scheme  for  the  Schrodinger  equation  we  will  discuss  an  accurate  discretization  of 
this  BC,  which  is  constructed  as  an  ABC  for  the  discrete  problem.  We  will  analyze  the  numerical  stability 
of  this  BC  and  present  numerical  tests.  An  extended  model,  formulated  on  the  level  of  density  matrices, 
allows  to  prescribe  incoming  wave  packets  at  the  contacts  in  combination  with  absorbing  the  outgoing 
waves.  In  the  classical  limit  (h  ->  0)  this  BC  corresponds  to  an  inflow  BC  for  a  kinetic  equation.  Finally 
we  will  discuss  ABC's  for  the  Wigner  equation.  Such  BC’s  have  been  used  for  the  numerical  simulation  of 
resonant  tunneling  diodes  via  the  relaxation-time  Wigner-Poisson  equation.  Most  of  these  studies  focused 
on  the  steady  state  behavior  of  the  device  for  an  applied  voltage.  The  known  ABC’s  for  the  Wigner 
equation,  however,  are  only  valid  for  short-time  transient  calculations  and  do  not  guarantee  the 
convergence  towards  the  correct  quantum  steady  state.  We  will  present  and  analyze  the  long-time 
modifications  of  these  ABC’s  for  the  relaxation-time  Wigner  equation  that  are  appropriate  for  the 
convergence  to  the  steady  state. 


WA2.  Modeling  Nonlinear  and  Chaotic  Dynamics  in  Semiconductor  Device  Structures,  Eckehard 
Scholl,  Institut  fur  Theoretische  Physik,  Technische  Universitdt  Berlin,  Hardenbergstr.  36,  10623  Berlin, 
Germany.  Semiconductors  represent  complex  nonlinear  dynamic  systems  whose  electrical  transport 
properties  may  exhibit  a  variety  of  instabilities  [1-3].  They  often  involve  switching  behavior,  self¬ 
generated  regular  or  chaotic  current  oscillations,  current  filamentation  and  solid-state  turbulence.  In  this 
paper  we  review  the  modeling  and  simulation  of  such  instabilities  with  a  special  emphasis  on  recent 
progress  in  the  application  to  semiconductor  microstructures.  Our  approach  involves  combinations  of 
semiclassical  Monte  Carlo  simulations  and  hydrodynamic  balance  equations.  In  the  center  of  interest  is  the 
nonlinear  spatio-temporal  dynamics  of  the  coupled  system  of  charge  carriers  and  electric  field.  The 
following  model  systems  are  treated  in  detail:  (i)  The  dynamics  of  current  filaments  in  the  regime  of  low- 
temperature  impurity  breakdown  in  p-Ge  and  n-GaAs  is  studied.  The  chaotic  behavior  of  laterally 
travelling  filaments  in  crossed  electric  and  magnetic  fields  is  considered  as  well  as  the  2D  simulation  of  the 
nascence  of  a  filament  upon  application  of  a  bias  voltage  [4].  The  transport  parameters  including  impact 
ionization  and  capture  rates  are  calculated  by  a  single  particle  Monte  Carlo  technique,  while  the  spatio- 
temporal  dynamics  is  obtained  self-consistently  from  the  carrier  continuity  and  Maxwell’s  equations,  (ii) 
Vertical  electrical  transport  in  layered  semiconductor  structures  like  the  heterostructure  hot  electron  diode 
(HHED)  is  considered  [5].  Periodic  as  well  as  chaotic  spatio-temporal  spiking  of  the  current  is  obtained 
and  analyzed  with  sophisticated  tools  from  nonlinear  dynamics.  In  particular  we  find  long  transients  of 
spatio-temporal  chaos  preceding  regular  spiking,  (iii)  Vertical  transport  in  a  superlattice  is  modeled  for 
high  fields  [6],  Self-generated  current  oscillations  or  the  formation  of  stable  stationary  electric  field 
domains  is  found.  We  demonstrate  that  this  behavior  is  strongly  affected  by  growth-related  imperfections 
and  weak  disorder. 
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WA3.  Studies  of  Chaotic  Transport  of  Electrons  in  Quantum  Boxes,  D.  K.  Ferry  and  G.  Edwards, 
Department  of  Electrical  Engineering,  Arizona  State  University,  Tempe,  AZ  85287-5706,  Recent  studies  of 
transport  through  ballistic  quantum  dot  resonators  have  revealed  a  complex  array  of  behavior,  including  the 
existence  of  ’’universal”  conductance  fluctuations  which  are  presumed  to  arise  not  from  disorder  scattering, 
but  from  the  chaotic  behavior  of  the  underlying  classical  dynamics  [1,2].  In  this  paper,  the  results  of  studies 
on  the  classical  ballistic  transport  of  carriers,  supposedly  within  a  quasi-two-dimensional  electron  gas, 
through  a  1.0  micron  square  structure  in  a  magnetic  field  are  presented.  The  boundary  on  one  side  is 
modified  in  a  manner  that  both  prohibits  straightthrough  trajectories.  The  presence  of  the  hard  walls, 
connecting  wires,  and  a  magnetic  field  are  thought  to  combine  to  introduce  chaos  in  the  structures.  Results 
for  B=0.01  T  are  shown  in  the  attached  figures.  We  will  also  discuss  the  effects  of  interelectronic  scattering 
in  the  fully  interacting  limit. 
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WA4.  Simulation  of  Quantum-Dot  Structures  in  Si/Si02,  Minhan  Chen  and  Wolfgang  Porod\ 
Department  of  Electrical  Engineering,  University  of  Notre  Dame,  Notre  Dame,  IN  46556.  We  present 
numerical  simulations  for  the  design  of  gated  few-electron  quantum  dot  structures  in  the  Si/Si02  material 
system.  The  motivation  for  this  work  is  to  investigate  the  feasibility  of  transferring  the  emerging 
technology  of  quantum  dot  fabrication  from  the  III-V  material  system,  where  it  was  pioneered  over  the  past 
few  years,  to  the  technologically  more  important  si/si02  structures.  Our  main  emphasis  is  on  the 
realization  of  coupled  quantum  dot  structures,  so-called  Quantum  Cellular  Automata  [1],  in  silicon.  Si 
appears  to  be  a  promising  candidate  due  to  the  excellent  insulating  behavior  of  thin  siOj  films  which  yields 
the  required  crisp  gate-control  of  the  potential  in  the  plane  of  the  two-dimensional  electron  gas  at  the 
Si/SiOj  interface.  In  our  simulations,  the  confining  potential  is  obtained  from  the  Poisson  equation  with  a 
Thomas-Fermi  charge  model.  Following  our  previous  work  [2,3],  the  electrostatic  potential  is  obtained 
both  in  the  semiconductor  domain  and  in  the  dielectric  above.  The  electronic  states  in  the  quantum  dot  are 
then  obtained  from  solutions  of  the  axially  symmetric  Schrodinger  equation.  We  explore  various  gate 
configurations  and  biasing  modes.  Our  simulations  show  that  the  number  of  electrons  can  be  effectively 
controlled  in  the  few  electron  regime,  especially  by  the  use  of  separately-biased  dual-gate  structures. 
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WAS.  Adiabatic  Switching  of  Quantum  Cellular  Automata,  Craig  5.  Lent  and  P.  Douglas  Tougaw, 
Department  of  Electrical  Engineering,  University  of  Notre  Dame,  Notre  Dame,  IN  46556.  We  have 
previously  proposed  a  way  of  using  coupled  quantum  dots  to  construct  digital  computing  elements.  Cells 
composed  of  a  few  quantum  dots  which  are  tunneling-coupled  are  occupied  by  few  electrons.  The 
Coulomb  interaction  within  the  cell  produces  two  possible  ground  states  with  quite  different  charge 
alignment  within  the  cell.  The  actual  ground  state  is  determined  by  the  Coulombic  (non-tunneling) 
interaction  with  neighboring  cells.  In  this  way  a  cellular-automata-like  architecture  can  be  built  up.  We 
have  shown  that  the  ground  state  of  an  inhomogeneous  cellular  array  can  be  mapped  onto  the  solution  of  a 
logical  or  arithmetic  problem  using  this  quantum  cellular  automata  (QCA)  approach  [1].  Here  we  will 
examine  the  problem  of  switching  such  QCA  arrays.  Abrupt  switching  has  been  discussed  in  the  past  and 
relies  on  intrinsic  dissipative  processes  to  relax  the  array  to  its  new  ground  state.  However,  a  clocked 
switching  mode  with  a  controlled  time-to-solution  may  be  desirable.  We  consider  here  switching  the  array 
slowly  so  that  it  is,  at  each  point  in  time,  in  the  instantaneous  ground  state.  We  will  present  results  of  direct 
time-dependent  solutions  of  the  few-electron  Schrodinger  equation  which  show  that  such  adiabatic 
switching  is  possible.  The  calculations  are  enabled  by  a  novel  use  of  operator-overloading  methods  in 
Fortran  90.  Implications  for  possible  future  device  architectures  using  adiabatic  QCA  arrays  are  discussed. 
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WA6.  Analysis  of  Q,-Independent  Single-Electron  Systems,  Alexander  N.  Korotkov  and  Konstantin  K 
Likharev,  Department  of  Physics,  State  University  of  New  York  at  Stony  Brook,  NY  11794-3800.  Fractional 
background  charge  Q^  is  a  major  problem  of  the  applied  single-electronics.  Most  of  single-electronic 
devices  suggested  so  far  (see,  e.g..  Ref.  1)  require  a  definite  value  of  Q^  with  the  margins  of  the  order  of 
±0.1e.  In  spite  of  evidence  of  relaxation  of  the  background  charge  in  some  experiments,  conditions  of 
such  relaxation  are  not  yet  clear.  The  goal  of  this  work  was  to  analyze  in  detail  two  single-electron  digital 
systems  which  may  operate  at  random  values  of  the  background  charge.  Both  systems  are  using  the 
periodic  dependence  of  the  current  in  the  single-electron  transistor  on  the  gate  voltage,  and  hence  its  ability 
to  perform  several  oscillations,  regardless  of  the  background  charge,  when  fed  by  a  gate  pulse.  These 
oscillations  may  be  easily  sensed  using  ordinary  FET  amplifier,  located  at  a  relatively  small  distance  L 
from  the  SET  transistor  (for  example,  for  a  realistic  output  SET  resistance  of  l00kQ>,a  ~  I  —  Gbl  S 
bandwidth  may  be  provided  2X  L<  100/l7t).  The  first  system  we  have  considered  is  a  DRAM,  were  a 
few-electron  charge  (--lOe)  may  be  written  into  the  initially  empty  trap  [1]  by  application  of  the  voltage 
pulse  to  the  trap  array.  Our  results  show  that  room  temperature  operation  and  density  in  excess  of  10" 
bits/cm^  may  be  achieved  in  such  RAMs,  using  -1-nm  nanolithography.  The  similar  principle  may  be 
applied  to  superdense  (up  to  -10*^  bit/cm^)  data  storage  in  a  layer  of  ultrafine  conducting  grains  (D~l  nm) 
separated  from  the  conducting  substrate  by  a  -10-nm-thick  insulator.  Binary  data  may  be  written  into  the 
grains  as  extra  few-electron  charges,  using  signal  voltage  applied  to  a  head  (tip)  moving  close  to  the 
surface.  Readout  of  the  data  may  be  performed  with  the  same  tip  carrying  the  SET/FET  transistor  pair 
discussed  above.  In  our  report,  methods  and  results  of  geometrical  analysis  and  numerical  simulation  of 
these  two  systems  will  be  presented  in  detail. 
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WA7.  Corrections  to  the  Capacitance  Between  Two  Electrodes  Due  to  the  Presence  of  a  Quantum 
Confined  System,  M.  Macucci,  Dipartimento  di  Ingegneria  dell'Informazione,  Via  Diotisalvi,  2  1-56126 
Pisa,  Italy;  and  Karl  Hess,  Beckman  Institute,  405  N  Mathews,  Urbana,  IL  61801.  For  structures 
enclosing  low-dimensional  systems  the  classical  concept  of  differential  capacitance  needs  to  be  augmented 
with  quantum  corrections  that  include  the  effects  of  confinement  energy.  Luryi  first  approached  the 
problem  by  deriving  the  quantum  capacitance  contribution  for  the  2D  electron  gas  of  an  infinite  MOS 
structure.  Our  aim  is  to  extend  this  concept  to  more  complex  situations,  characterized  by  additional 
degrees  of  confinement,  and  to  determine  whether  quantum  corrections  may  play  a  significant  role  in 
scaled-down  classical  devices.  Our  approach  is  based  on  the  evaluation  of  the  free  energy  of  the  combined 
device  plus  voltage  source  system,  a  method  often  used  in  Coulomb  Blockade  studies.  By  minimizing  the 
free  energy,  we  can  determine  the  number  of  electrons  in  the  quantum  system  (2DEG,  quantum  dot,  etc.) 
for  a  given  applied  voltage  and  the  thresholds  for  the  addition  of  each  electron.  A  few  simple  cases  can  be 
treated  analytically,  and  we  show,  as  an  example,  that  Luryi's  formula  can  be  retrieved.  The  case  of  the  2D 
electron  gas  is  particularly  simple,  because  of  the  quadratic  dependence  of  the  confinement  energy  on  the 
number  of  electrons:  this  makes  possible  to  consider  an  equivalent  circuit  including  a  "quantum  capacitor”. 
If  further  confinement  is  added  such  as  in  the  case  of  a  quantum  dot,  this  quadratic  dependence  disappears, 
and  it  is  not  possible  to  resort  to  simple  equivalent  circuits.  A  numerical  approach,  which  also  allows  the 
study  of  more  realistic  models  becomes  necessary.  We  solve  the  Schrodinger  and  Poisson  equations  self- 
consistently  in  structures  with  cylindrical  symmetry,  thus  obtaining  the  electron  eigenvalues  and 
eigenfunctions  in  the  quantum  dot  and  the  charge  induced  on  the  electrodes,  from  which  the  free  energy 
and  thereby  the  tunneling  thresholds  of  each  electron  can  be  computed.  The  differential  capacitance  is 
finally  derived  by  taking  the  ratio  of  the  variation  of  charge  induced  on  the  electrodes  to  that  of  the  bias 
voltage. 


WAS.  Simulation  of  Dissipative  Quantum  Transport  in  Electron  Distributed  Bragg  Reflectors, 
Leonard  F.  Register  and  Karl  Hess,  Beckman  Institute  and  Coordinated  Science  Laboratory,  University  of 
Illinois  at  Urbana-Champaign,  61801.  The  effects  of  inelastic  scattering  on  carrier  transport  through 
electron  distributed  Bragg  reflectors  (electron  DBRs)'’^  is  simulated  using  both  "Schrodinger  equation 
(based)  Monte  Carlo"  (SEMC)”  and,  more  simply,  a  complex  absorbing  potential  approach  to  scattering 
[5].  Analogous  to  optical  DBRs,  electron  DBRs  employ  interference  among  the  relatively  weak  partial 
reflections  of  the  electron  wave  function  at  each  heterointerface  to  achieve  an  overall  high  reflectivity. 
However,  inelastic  scattering  processes  in  semiconductors  limit  the  phase  coherence  length  of  the  electron 
wave  function,  and,  thus,  can  be  expected  to  limit  the  reflectivity  of  electron  DBRs.  SEMC  and  the 
complex  "absorbing"  potential  model  allow  simulation  of  such  coherence  limiting  dissipative  quantum 
transport.  SEMC  is  specifically  formulated  to  allow  first-principles  simulation  of  dissipative  quantum 
transport  and  is  rigorously  quantum  mechanical.  However,  the  numerical  algorithm  has  much  in  common 
with  semiclassical  Monte  Carlo  methods.  For  carrier-phonon  interactions,  a  set  of  Schrodinger  equations  is 
solved  simultaneously  (and  deterministically)  for  the  carrier,  with  the  individual  equations  corresponding 
to  the  discrete  initial  or  "trunk"  state  and  various  final  or  "branch"  states  of  the  phonon  system,  with 
energies  separated  by  plus  or  minus  the  phonon  energy  for  phonon  absorption  and  emission,  respectively. 
Only  the  non-local  coupling  potentials  between  the  trunk  and  branch  states  are  obtained  by  Monte  Carlo 
sampling,  of  (the  spatial  correlation  functions  of)  the  carrier-phonon  interactions  [3].  In  contrast,  the  use  of 
absorbing  potentials  within  Schrodinger’s  equation  in  lieu  of  the  branch  states  to  represent  scattering, 
allows  more  simple  and  efficient  if  less  rigorous  modeling  of  the  loss  of  phase  coherence  [5]. 
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WPAl.  The  Coupled  Optoelectronic  Problems  of  Quantum  Well  Laser  Operation,  Karl  Hess  and  Matt 
Grupen,  The  Beckman  Institute,  Department  of  Electrical  and  Computer  Engineering,  University  of 
Illinois,  405  N.  Mathews  Ave.,  Urbana,  IL  6I80I.  The  optical  output  of  a  quantum  well  laser  depends  on 
several  different  physical  processes.  These  include  classical  carrier  transport,  mesoscopic  transport  over  a 
quantum  structure,  intrasubband  scattering,  quantum  confinement,  radiative  recombination,  and 
electromagnetic  wave  guiding.  Each  process  is  important  to  laser  operation,  and  each  poses  a  complex 
theoretical  problem.  The  different  optics  and  electronics  related  problems  associated  with  laser  operation 
have  been  treated  separately  and  independently  in  great  detail.  Simulators  have  been  developed  to  calculate 
classical  transport  throughout  an  entire  two-dimensional  laser  cross  section  [1].  Other  sophisticated  codes 
have  been  developed  to  determine  the  reflection  and  transmission  of  carriers  incident  on  a  quantum  well 

[2]  and  the  scattering  dynamics  of  bound  carriers  [3].  Still,  other  research  has  been  conducted  on  the  band 
structure  within  a  quantum  well  and  the  matrix  elements  associated  with  optical  transitions  [4].  Also, 
considerable  work  has  been  done  on  the  accurate  solution  of  Maxwell's  equations  in  a  dielectric  wave 
guide  [5]  and  on  emission  from  a  microcavity  [6].  However,  these  separate  treatments  cannot  reveal  how 
these  different  problems  couple  and  how  their  interaction  affects  laser  operation.  A  laser  simulator 
Minilase  II  has  been  developed  to  address  all  of  the  aforementioned  problems  self-consistently  and 
completely  coupled  within  the  limits  of  workstation  computational  capacity.  This  paper  will  describe  the 
ways  in  which  these  problems  are  modeled  and  the  methods  used  to  couple  them.  It  will  also  present 
certain  results  that  highlight  coupling  between  different  physical  processes  and  which,  therefore,  can  only 
be  obtained  through  such  a  self-consistent  approach. 
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WPA2.  Modeling  Light-Extraction  Characteristics  of  Packaged  Light-Emitting  Diodes,  D.  Z.-Y.  Ting 
and  T.  C.  McGill,  Thomas  J.  Watson,  Sr.  Laboratory  of  Applied  Physics,  Mail  Stop  128-95,  California 
Institute  of  Technology,  Pasadena,  California  91125.  We  employ  a  Monte  Carlo  ray  tracing  technique  to 
model  light-extraction  characteristics  of  light-emitting  diodes.  By  relaxing  restrictive  assumptions  on 
photon  traversal  history,  our  method  improves  upon  available  analytical  models  for  estimating  light- 
extraction  efficiencies  from  bare  LED  chips,  and  enhances  modeling  capabilities  by  realistically  treating 
the  various  processes  which  photons  can  encounter  in  a  packaged  LED.  Our  method  is  not  only  capable  of 
calculating  extraction  efficiencies,  but  can  also  provide  extensive  statistical  information  on  photon 
extraction  processes,  and  predict  LED  spatial  emission  characteristics.  Preliminary  results  from  a  sample 
parametric  study  indicate  that  simulations  using  our  method  can  be  performed  very  rapidly  on  modem 
workstations.  We  expect  that,  with  some  refinement,  our  method  could  be  used  as  a  computer-aided  design 
tool  for  optimizing  light-extraction  efficiencies,  and  for  designing  LED  spatial  emission  characteristics. 


WPA3.  Gain  Calculation  in  a  Quantum  Well  Laser  Simulator  Using  an  Eight  Band  k  .  p  Model,  F. 

Oyafuso,  P.  Von  Allmen,  M.  Grupen,  K.  Hess,  Beckman  Institute,  University  of  Illinois,  Urbana,  IL  61801. 
Gain  spectra  for  a  strained-layer  GagIn^As/Al.iGa9As  quantum  well  laser  are  calculated  using  an  eight 
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band  k  .  p  model.  Effects  of  non-parabolicity  of  the  energy  dispersion  are  entered  in  a  QW  laser  simulator 
(MINILASE-II).  The  k .  p  calculation  yields  the  bandstructure  and  momentum  matrix  elements  for  a 
superlattice  grown  along  the  [001]  direction.  The  calculation  works  within  the  envelope  function 
approximation,  where  the  envelope  functions  are  expanded  in  a  plane  wave  basis.  The  barrier  between  the 
wells  of  the  superlattice  is  made  sufficiently  large  so  that  the  bound  state  solutions  are  expected  to  match 
those  of  a  single  quantum  well.  The  k .  p  results  are  then  exported  to  MINILASE-II  in  the  form  of  a 
density  of  states  and  an  energy-dependent  averaged  momentum  matrix  element.  The  matrix  elements  are 
averaged  over  the  entire  ID  mini-band  Brillouin  zone  and  the  set  of  in-plane  wave  vectors  corresponding  to 
a  given  pair  of  initial  and  final  energies.  For  each  transition  energy,  many  pairs  of  conduction  and  valence 
band  states  can  contribute  because  of  the  anisotropy  and  multivalued  nature  of  the  band  structure.  A 
broadening  factor  is  introduced  within  the  laser  simulator.  To  reduce  computation  time,  the  calculations  are 
coupled  to  Poisson's  equation  within  MINILASE-II.  In  other  words,  a  fiat  band  approximation  is  assumed 
for  the  k .  p  calculation,  and  the  resulting  density  of  states  is  pinned  to  the  lowest  subbands  calculated 
within  MINILASE-II.  Results  will  be  presented  for  the  gain  and  modulation  response  using  this  model  for  a 
TE  optical  field  and  will  be  compared  with  the  results  obtained  from  a  simple  parabolic  Kane  model  with  a 
constant  matrix  element. 


WPBl.  Moving  Adaptive  Unstructured  3-D  Meshes  in  Semiconductor  Process  Modeling 
Applications,  Harold  Trease,  Andrew  Kuprat,  Denise  George,  Eldon  Linnebur,  Los  Alamos  National 
Laboratory,  and  R.  Kent  Smith,  AT -¥7  Bell  Laboratories.  The  next  generation  of  semiconductor  process 
and  device  modeling  codes  will  require  3-D  mesh  capabilities  including  moving  volume  and  surface  grids, 
adaptive  mesh  refinement  and  adaptive  mesh  smoothing.  To  illustrate  the  value  of  these  techniques,  a  time 
dependent  process  simulation  model  was  constructed  using  analytic  functions  to  return  time  dependent 
dopant  concentration  and  time  dependent  SiO^  volume  and  surface  velocities  .  Adaptive  mesh  refinement 
and  adaptive  mesh  smoothing  techniques  were  used  to  resolve  the  moving  boron  dopant  diffusion  front  in 
the  Si  substrate.  The  lower  SiO^  surface  was  modeled  as  a  moving  surface  grid  which  advanced  through 
the  background  grid,  converting  the  polysilicon  and  silicon  in  its  path  to  SiOj.  An  interesting  problem 
addressed  in  the  moving  surface  implementation  of  the  oxide  expansion  was  maintaining  material 
interfaces  by  reapplying  the  material  region  definitions  at  each  time  step.  The  expanding  volume  of  SiOj, 
in  turn,  pushed  the  unreactive  nitride  mask  ahead,  requiring  that  the  nitride  region  be  treated  as  an 
incompressible  volume.  The  novel  adaptive  mesh  smoothing  technique  for  resolving  the  evolving  dopant 
concentration  in  the  underlying  Si  substrate  involves  minimizing  an  error  functional  over  the  volume  mesh. 
The  error  functional  chosen  is  the  L^  norm  of  the  gradient  of  the  error  between  the  true  dopant 
concentration  and  the  piecewise  linear  approximation  over  the  tetrahedral  mesh.  By  minimizing  the 
gradient  of  the  error  (rather  than  just  the  error)  we  assure  that  the  mesh  is  optimal  for  representing  evolving 
solution  gradients.  This  implies  correct  mesh  alignment  with  the  dopant  field  as  well  as  high  quality 
tetrahedral  shapes.  Also  implemented  is  constrained  boundary  smoothing,  wherein  the  moving  SiO/Si 
interface  is  represented  by  moving  nodes  that  correctly  track  the  interface  motion,  and  which  use  their 
remaining  degrees  of  freedom  to  minimize  the  aforementioned  error  functional.  This  implies  optimal 
tetrahedral  shape  and  alignment  is  obtained  even  in  the  neighborhood  of  a  moving  boundary.  At  each  time 
step,  a  topological  “reconnection”  step  assures  that  no  negative  coupling  coefficients  are  present  in  the 
mesh.  The  combination  of  adaptive  refinement,  adaptive  smoothing,  and  mesh  reconnection  give  excellent 
front  tracking,  feature  resolution,  and  grid  quality  for  finite  volume/finite  element  computation. 


WPB2.  Computational  Aspects  of  Topography  Evolution  During  Low  Pressure  Deposition  and  Etch 
Processes,  Timothy  S.  Cale,  Center  for  Solid  State  Electronics  Research,  Department  of  Chemical,  Bio 
Materials  Engineering,  Arizona  State  University,  Tempe,  AZ  85287^6206.  Low  pressure  deposition  and 
etch  steps  are  common  in  device  fabrication.  Transport  and  reaction  at  low  pressures  in  micron  scale 
features,  where  species  transport  through  the  vacuum  is  essentially  collisionless,  is  well  modeled  by  the 
'ballistic  transport  and  reaction  model'  (BTRM)  [1,2].  The  integrodifferential  equations  which  represent 
the  BTRM  include  terms  which  express  how  the  ballistic  species  interact  with  the  surface;  i.e.,  adsorb  onto, 
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react  with  or  re-emit  from  the  surface.  The  solution  of  the  BTRM  equations  yields  the  distributions  of  the 
fluxes  of  ballistic  species  and  surface  coverages  by  physisorbed  and/or  chemisorbed  species  over  the 
surface.  In  turn,  these  yield  rates  of  surface  evolution.  Surface  evolution  is  performed  using  the  method  of 
characteristics  [3].  In  order  to  perform  accurate,  predictive  topography  simulation,  submodels  for  the 
reaction  kinetics,  ballistic  species  flux  distributions  and  surface  diffusion  of  adsorbed  species  must  be 
established.  Although  simulations  of  the  processing  of  most  device  structures  require  complete  three 
dimensional  transport  and  surface  movement  algorithms  [4],  we  have  focused  on  structures  for  which  three 
dimensional  transport  and  two  dimensional  surface  profile  evolution  is  appropriate.  Infinitely  long 
trenches  or  lines  and  features  of  circular  horizontal  cross  section  (idealized  vias  or  contact  holes)  fall  into 
this  category.  These  structures  are  ideal  to  establish  models  for  species  transport  and  chemical  reactions, 
because  film  profiles  can  be  experimentally  determined  using  electron  microscopy.  Models  can  be 
developed  and  validated  by  comparing  simulated  and  experimental  film  profiles.  These  models  can  then 
be  used  with  confidence  in  deposition  and  etch  process  simulations  for  fully  three  dimensional  surfaces.  In 
this  paper,  I  will  review  1)  the  BTRM,  2)  the  methods  used  to  solve  the  governing  equations  of  the  BTRM, 
with  and  without  surface  diffusion,  and  3)  the  implementation  of  the  method  of  characteristics  for  multiple 
material  substrates. 


[1]  T.S.  Cale  and  G.B.  Raupp,  J.  Vac.  Sci.  Tech.  B8(4),  649(1990). 

[2]  T.S.  Cale  and  G.B.  Raupp,  J.  Vac.  Sci.  Tech.  B8(6),  1242(1990). 

[3]  D.  Ross,  J.  Electrochem.  Soc.  135(5)  1235  and  1260  (1988). 

[4]  H.  Liao  and  T.S.  Cale,  Thin  Solid  Films  236(1),  352  (1993). 


WPB3.  Hierarchical  Process  Simulation  for  Nanoelectronics,  Robert  W.  Dutton,  Stanford  University 
Center  for  Integrated  Systems,  Stanford,  CA  94305.  Over  a  period  of  more  than  two  decades  the  field  of 
process  simulation  for  integrated  circuits  has  become  established  as  an  essential  enabling  technology.  Over 
this  same  period  the  critical  dimensions  for  devices  have  steadily  moved  from  the  regime  of  10  pm  to  well 
below  1.0  pm  -  the  most  recent  technology  papers  are  pushing  the  limits  of  0.1  pm  devices.  With  dielectric 
layers  below  10  nanometers  and  junction  depths  below  100  nanometers,  it  is  certainly  safe  to  say  that  we 
are  now  entering  the  era  of  nanoelectronics.  The  focus  of  this  talk  will  be  to  look  critically  at  the  recent 
history  of  developments  in  process  simulation  and  to  consider  the  opportunities  and  challenges  facing  this 
new  era  of  IC  technology  development.  The  following  gives  an  outline  of  several  topics  and  results  to  be 
discussed.  Bulk  process  simulation  relates  to  transformations  of  the  silicon  wafer  to  create  the  active 
devices  and  it  continues  to  pose  many  challenges  both  for  simulation  and  characterization.  The  topics 
related  to  ion  implantation,  thermal  diffusion  and  oxidation  will  each  be  considered  both  in  terms  of  atomic 
level  effects  as  well  as  based  on  the  more  traditional  continuum  models.  There  is  a  growing  need  for  and 
emphasis  on  first-principles  modeling  of  these  bulk  processing  effects,  exactly  because  of  the  nanometer 
scale  effects  and  device  requirements.  Even  where  continuum  modeling  can  be  applied,  there  are  needs  to 
extract  model  parameters  and  understand  the  interdependence  In  one  sense  these  issues  very  much  parallel 
those  faced  in  developing  the  hierarchy  of  device  simulation  models;  both  the  similarities  and  differences 
will  be  considered.  Simulation  of  multi-level  interconnections,  the  back-end  of  the  process,  has  become  a 
major  topic  of  focus  and  concern.  Both  the  cost  and  performance  of  technology  is  critically  impacted  by 
interconnect  technology.  In  contrast  to  the  front-end  or  bulk  processing,  the  materials  and  process  physics 
put  a  very  different  set  of  requirements  on  the  simulators— greater  emphasis  on  surfaces  and  kinetics  of 
deposition  and  etching.  While  the  development  of  such  back-end  simulators  is  many  years  behind  that  of 
front-end  simulators,  there  are  new  opportunities  to  share  both  atomic-scale  modeling  techniques  and 
software  infrastructure.  Examples  will  be  presented  that  show  very  exciting  progress  to  fully  bridge  the 
hierarchy  from  atomic  through  continuum  modeling  techniques  in  this  area.  Finally,  the  software 
engineering  and  other  computational  issues  related  to  hierarchical  process  simulation  will  be  considered. 
Deep  software  reconfigurability  in  support  of  new  and  changing  problems  and  numerical  techniques  is 
needed  to  ensure  both  integrity  and  timely  development  on  software.  With  the  growing  computational 
demands  arising  from  complex  physical  definition  and  full  3D  treatment,  software  in  the  simulation 
hierarchy  needs  to  exploit  possible  parallel  and  scalable  techniques  including  solver  libraries  and 
geometry/grid  servers.  With  our  experience  with  parallelization  of  the  device  simulation  code,  it  can  be 
projected  that  effective  use  of  computers  with  hundreds  of  processors  will  achieve  performance  levels  that 
can  support  engineering  application  of  a  full  hierarchy  of  process  modeling  capabilities. 
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WPB4.  A  Compound  Semiconductor  Process  Simulator  and  its  Application  to  Mask  Dependent 
Undercut  Etching,  Masami  Kumagai  and  Kiyoyuki  Yokoyama,  NTT  Opto-electronics  Laboratories,  and 
Satoshi  Tazawa,  NTT  LSI  Laboratories.  A  process  simulator  for  the  compound  semiconductor  process  is 
developed.  This  simulator  is  designed  for  describing  the  etching  and  deposition  processes  of  compound 
semiconductors.  A  block  diagram  of  the  simulator  is  given  in  Figure  1.  It  consists  of  three  user- interface 
modules  and  three  calculation  modules.  The  calculation  modules,  the  "Gas  Distribution  Module",  the 
"Crystal  Orientation  Module"  and  the  "Surface  Reaction  Module",  determine  the  time  evolution  of  the 
etching  and  deposition  processes.  The  "Gas  Distribution  Module"  and  the  "Crystal  Orientation  Module" 
supplement  the  "Surface  Reaction  Module"  and  enhance  its  applications.  In  compound  semiconductor 
processes,  especially  in  optical  device  processes,  crystal  orientation  dependence  plays  an  important  role. 
Compound  semiconductors,  such  as  GaAs  or  InP,  have  anions  and  cations  while  Si  is  made  of  a  single 
atomic  species.  This  binary  nature  yields  a  variety  of  processed  profiles  including  mesa  or  reverse  mesa 
shape.  In  the  optical  device  processes,  the  mesa  and  the  reverse  mesa  shape  formation  is  very  common  and 
thus  detailed  representations  are  required  for  the  process  simulation.  Figures  2a  and  2b  give  the  simulated 
etching  profiles  for  mesa  and  reverse  mesa  faces,  respectively.  Mesa  and  reverse  mesa  shape  formation  is 
quite  sensitive  to  process  conditions.  Figures  3a  and  3b  show  the  experimental  wet  etching  profiles  for  an 
InP  substrate  with  two  different  mask  materials,  Si02  and  SiN,  respectively.  The  Si02  masked  InP  shows 
the  large  undercut  etching  and  the  side  faces  become  throughly  mesa  shaped,  while  the  SiN  masked  InP 
gives  lesser  undercut  etching  and  the  mesa  shaped  portions  are  significantly  shrunk.  This  etching  process 
was  simulated  by  using  the  virtual  metamorphic  layer  model,  which  employs  a  thin  virtual  layer  between 
the  mask  material  and  the  InP  substrate  and  assumes  this  layer  is  changed  during  the  mask  deposition 
process.  The  simulation  results  are  presented  in  Figs.  4a  and  4b.  Good  agreement  between  the  Fig.3  and 
Fig.  4  can  be  seen.  The  simulation  conditions  for  4a  and  4b  are  equivalent  except  for  the  etching  rate  of 
each  metamorphic  layer.  The  etching  rate  of  the  metamorphic  layer  under  the  Si02  mask  (Fig.  4a)  is 
assumed  to  be  five  times  that  of  the  substrate  material,  while  that  under  the  SiN  mask  (Fig.  4b)  is  assumed 
to  be  the  same  as  the  substrate  material.  The  mask  etching  rate  is  taken  to  be  zero.  In  summary,  the 
process  simulator  describes  the  complicated  etching  processes  on  compound  semiconductors  yielding 
information  about  the  physical  properties  of  the  processes  as  well  as  the  process  conditions  for  the  desired 
processed  shapes.  Details  of  the  simulation  procedures  and  further  results  will  be  presented  by  an  oral 
presentation. 


WPB5.  The  Combination  of  Equipment  Scale  and  Feature  Scale  Models  for  Chemical  Vapor 
Deposition  via  a  Homogenization  Technique,  M.  Gobbert*,  T.S.  Cale#,  and  C.  Ringhofer*,  *Department 
of  Mathematics,  Arizona  State  University,  Tempe,  AZ  85287-1804;  #Center  for  Solid  State  Electronics 
Research,  Dept,  of  Chemical,  Bio  &  Materials  Engineering,  Arizona  State  University,  Tempe,  AZ  85287- 
6206.  The  conformality  and  composition  of  deposited  films  on  the  scale  of  device  features  impact  device 
performance  and  reliability.  Variations  of  these  properties  along  the  wafer  during  processing  lead  to 
unwanted  variations  in  device  performance.  Thus,  an  acceptable  deposition  process  provides  both  average 
values  and  across  wafer  uniformity  of  conformality  and  composition  within  process  dependent,  specified 
ranges  in  addition  to  an  acceptable  wafer  throughout.  In  recent  years,  equipment  scale  models  have 
successfully  been  used  to  establish  acceptable  operating  conditions  for  existing  reactors  and  to  design 
reactors  which  provide  the  necessary  across  wafer  uniformity  [1].  More  generally,  feature  scale 
simulations  should  be  coupled  with  equipment  scale  simulations  in  order  to  ensure  acceptable  response, 
i.e.,  to  predict  film  properties  and  uniformity  as  functions  of  reactor  setpoints.  In  previous  studies,  reactor 
scale  simulators  have  teen  used  to  provide  the  boundary  conditions  (fluxes)  for  feature  scale  simulators 
[2].  This  approach  assumes  that  feature  scale  processes  do  not  influence  the  local  species  concentrations  at 
the  wafer  surface;  i.e.,  increases  in  the  local  deposition  area  due  to  topography  do  not  affect  film  properties. 
If  this  is  not  the  case,  then  the  predictions  of  local  film  properties  can  be  improved  upon  by  using 
information  from  the  feature  scale  simulator  in  the  reactor  scale  simulator.  A  truly  global  simulator  would 
solve  feature  scale  and  reactor  scale  models  simultaneously  and  account  for  local  'loading'  effects  due  to 
the  presence  of  topography,  when  appropriate.  In  this  paper,  we  present  a  method  for  combining  the  vastly 
different  length  scales  appropriate  for  equipment  and  feature  scale  simulators  based  on  a  homogenization 
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technique  from  asymptotic  analysis.  The  reactor  scale  model  is  coupled  to  the  surface  via  a  near  surface 
region.  The  near  surface  region  accounts  for  surface  topography  using  asymptotic  analysis  to  account  for 
the  micro  structure  of  the  wafer  surface.  The  equations  which  govern  transport  and  deposition  inside  the 
features  are  solved  simultaneously,  and  the  appropriate  species  fluxes  will  be  included  in  the  boundary 
conditions  of  the  near  surface  region,  which  are  in  turn  will  be  used  as  boundary  conditions  for  the 
equipment  scale  equations.  This  model  is  capable  of  predicting  the  effects  of  feature  aspect  ratios  and 
feature  spacing  on  film  properties.  This  approach  has  been  validated  mathematically  and  numerically  [3]. 
We  explain  the  model  and  discuss  predictions  of  our  simulator  based  on  the  model  by  analyzing  the 
deposition  of  silicon  dioxide  from  teu-aethoxysilane  (TEOS).  In  particular,  we  will  discuss  the  effects  of 
pattern  density  variations  on  local  film  thickness  uniformity. 


[1]  D.W.  Studiner,  J.T.  Hillman,  R.  Aurora  and  R.F.  Foster,  in  Advanced  Metallization  for  ULSI 
Applications  1992,  T.  S.  Cale  and  F.  S.  Pintchovski,  eds..  Materials  Research  Society,  1993.  211. 

[2]  T.S.  Cale,  J.-H.  Park,  T.H.  Gandy,  G.B.  Raupp  and  M.K.  Jain,  Chemical  Engineering  Communications, 
119,  197(1993). 

[3]  M.K.  Gobbert  and  C.A.  Ringhofer,  An  asymptotic  analysis  for  a  model  of  chemical  vapor  deposition  on 
a  micro  structured  surface,  submitted. 
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